32.sagews

# profit on a dominant assurance contract
K,V = var("K,V")
N = 100
sigma, mean = 1/12, 0.5

def G(x):
    return (1+erf((x-mean)/(sigma*sqrt(2.0))))/2

G(0.2)

0.000159108590157553

def z(V,K):
    return V*K*binomial(N,floor(K))*(1-G(V))^K*G(V)^(N-K)

z(0.5,40)

0.216877334232760

parametric_plot3d([K, V, z], (K,0, 100), (V, 0, 1), plot_points=[101,51], mesh=True)