# Solving a numeric statistics problem - R [closed]

I am quite stuck with solving some complicated numerical equation

I would like to solve the following equation:

$(1-k)\tilde{\alpha}+kf(\tilde{\alpha})=C$

where $0<k<1$ and $f(\tilde{\alpha})$ is a right tail percentage of a multivariate distribution of a length n with the covariance matrix is 1 on the diagonal and n-2/n-1 in all the other entrences and the mu vector equals zero in all n entrances.

then:

$f(\tilde{\alpha})=P(Z_1>z_{1-\tilde{\alpha}},...,Z_n>z_{1-\tilde{\alpha}})$

when $n=\infty$ we get that $\tilde\alpha=C$

generally:

$\tilde{\alpha}^n\leq f(\tilde{\alpha})\leq\tilde{\alpha}$

and $C\leq\tilde{\alpha}\leq1$

Thanks a lot.

• Is it a multivariate normal distribution or something else? Please confirm if $\tilde{\alpha}$ is the unknown? What does $m$ and $\alpha$ represent? Are they simply constant given to you? Do you know how to evaluate $f(\tilde{\alpha})$ or is this the main difficulty you are faced with? May 15 '16 at 18:32

From what I understand, you are basically interested in solving the equation numerically, over all possible $\tilde{\alpha}$. In case you are OK with Python:

from scipy.optimize import fmin_l_bfgs_b as bfgs
from scipy.stats import norm
from scipy.stats import mvn
import numpy as np
k = 0.2
c = 0.05
n = 3

def get_z(alpha):
return norm.ppf(alpha)

def get_sigma():
sigma = np.ones((n,n))*(n-2)/(n-1)
return sigma+np.identity(n)*1.0/(n-1)

def get_f(alpha):
mu = np.zeros(n)
S = get_sigma()
low = np.ones(n)*get_z(1-alpha)
upp = np.ones(n)*1000
p, i = mvn.mvnun(low, upp, mu, S)
return p

def g(alpha):
return (k*alpha + (1-k)*get_f(alpha)-c)**2

def optimize():