<jupyter_start><jupyter_text>1) von mises packages   - dvonmises<jupyter_code>def d_vonmises(x, mu, kappa, log = False):
    if (type(x) != list):
        x=[x]
    pdf = np.zeros(len(x))
    if (log):
        if (kappa < 100000):
            pdf = -(np.log(2*math.pi)+np.log(scp.ive(0, kappa)) + kappa) + kappa*(np.cos(np.subtract(x - mu)))
        else:
            if (((x-mu)%(2*math.pi))==0):
                pdf = math.inf
            else:
                pdf = -math.inf
    else:
        if (kappa == 0):
            pdf = np.repeat(1/(2*np.pi), len(x))
        elif (kappa < 100000):
            pdf = 1/(2 * np.pi * scp.ive(0, kappa)) * (np.exp(np.subtract(np.cos(np.subtract(x, mu)), 1)))**kappa
        else:
            if (np.mod(np.subtract(x, mu),(2*np.pi))==0):
                pdf = math.inf
            else:
                pdfm = 0
    return(pdf)
d_vonmises(1,1,6)
%timeit d_vonmises(1,1,6)<jupyter_output>22.4 µs ± 607 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)<jupyter_text>- qvonmises<jupyter_code>def q_vonmises(p, mu = 0 ,  kappa = None, from_ = None, tol = np.finfo(float).eps**0.6):
       
    if (type(p) != list):
        p = np.array([p])
    else:
        p = np.array(p)
        
    epsilon = 10 * np.finfo(float).eps  
    
    if (np.any(p > 1) or np.any(p < 0)): 
        raise ValueError("p must be between [0,1]")
        
    if (pd.isnull(from_)): 
        from_ = mu - np.pi
   
    n = len(p)
    mu = (mu - from_)%(2 * np.pi)      ## from is a keyword    
    
    if (len([mu]) != 1): 
        raise ValueError("is implemented only for scalar mean")   
        
    if (pd.isnull(kappa)): 
        raise ValueError("kappa should be provided")   
        
        
    def zeroPvonmisesRad(x, p, mu, kappa):
        if (np.isnan(x)):         
            y = np.nan              
        else: 
            integration = integrate.quad(lambda x: d_vonmises(x, mu, kappa), 0, x)
            y = integration[0] - p         ##integration[0] will give the value
        return(y);
    
    
    value = np.repeat(np.nan, p.size)
    for i in range(p.size):
        try:
            value[i] = root(lambda x: zeroPvonmisesRad(x, p[i], mu, kappa), 0, 2 * np.pi - epsilon)
        except:
            pass
            if(p[i] < (10 * epsilon)):
                value[i] = 0
            elif (p[i] > (1 - 10 * epsilon)):
                value[i] = 2 * np.pi - epsilon         
    value += from_
    return(value)

q_vonmises(0.5,1,6)
%timeit q_vonmises(0.5,1,6)<jupyter_output>6.36 ms ± 329 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)<jupyter_text>- pvonmises<jupyter_code>def p_vonmises(q, mu, kappa, tol = 1e-020):
    from_ = mu - np.pi
    mu = (mu - from_) % (2 * np.pi)
    
    
    if (type(q) != list):
        q = [q]
    q = np.mod(np.subtract(q, from_), (2 * np.pi))
    q = np.mod(q,(2 * np.pi))
    n = len(q)
    mu = mu % (2 * np.pi)
    def fn_mu0(q,kappa,tol):
        flag = 1
        p = 1
        sum_ = 0
        while(flag):
            term = (iv(p, kappa) * np.sin(np.multiply(q, p)))/p
            sum_ = sum_ + term
            p = p + 1
            if (abs(term) < tol):
                flag = 0
        return(np.divide(q,(2 * np.pi)) + sum_/(np.pi * iv(0, kappa)))
    result = np.repeat(np.nan, n)
    if (mu == 0):
        for i in range(0,n):
            result[i] = fn_mu0(q[i], kappa, tol)
    else:
        for i in range(0,n):
            if (q[i] <= mu):
                upper = (q[i] - mu) % (2 * np.pi)
                if (upper == 0):
                    upper = 2 * np.pi
                lower = (-mu) % (2 * np.pi)
                result[i] = fn_mu0(upper, kappa, tol) - fn_mu0(lower, kappa, tol)
            else:
                upper = q[i] - mu
                lower = mu % (2 * np.pi)
                result[i] = fn_mu0(upper, kappa, tol) + fn_mu0(lower, kappa, tol)
    return(result)
p_vonmises([2,0.8],2,6)
%timeit p_vonmises([2,0.8],2,6)<jupyter_output>811 µs ± 38.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)<jupyter_text>- rvonmises<jupyter_code>def r_vonmises(n,mu,kappa):
    a = np.random.uniform(0,1,n)
    b = [q_vonmises(x,mu,kappa) for x in a]
    c = np.squeeze([l.tolist() for l in b])
    return(np.array([a%(2*np.pi) for a in c]))
rv_ls = r_vonmises(50,10,3)
n=50
rv_ls
a = (np.cos(rv_ls)+ np.random.normal(scale=0.05,size=n))
b = (np.sin(rv_ls)+ np.random.normal(scale=0.05,size=n))
plt.plot(a, b,'o', color='black',alpha=0.3)<jupyter_output><empty_output><jupyter_text>2) Mixture Distributions of von mises Fisher  a) Method 1 - Only p=2 von mises distribution (input - angles)<jupyter_code>a=list(range(1,11))
b=list(range(11,21))
dat_ls=[math.atan(b/a) for a,b in (list(zip(a,b)))]

c=[10,15,20,25,100,110,120,130,140]

dat_ls = [x*pi/180 for x in c]
dat_ls
dat_vector = [[math.cos(x),math.sin(x)] for x in dat_ls]
p=len(dat_vector[0])
dat_vector

knum = 2
n = len(dat_ls)
mu=[1,0.5]
kappa=[10,10]
p_i = [0.3,0.7] 
def vmf_pdf(x, mu, kappa, log = False):
    if (type(x) == int):
        x = [x]
    if (type(x) == float):
        x = [x]
    vm = np.zeros(len(x))
    if (log):
        if (kappa == 0):
            vm = np.log(np.repreat(1/(2*pi), len(x)))
        elif (kappa < 100000):
            vm = -(np.log(2*math.pi)+np.log(scp.ive(0, kappa)) + kappa) + kappa*(np.cos(np.subtract(x - mu)))
        else:
            if (((x-mu)%(2*math.pi))==0):
                vm = math.inf
            else:
                vm = -math.inf
    else:
        if (kappa == 0):
            vm = np.repeat(1/(2*np.pi), len(x))
        elif (kappa < 100000):
            vm = 1/(2 * np.pi * scp.ive(0, kappa)) * (np.exp(np.subtract(np.cos(np.subtract(x, mu)), 1)))**kappa
        else:
            if (np.mod(np.subtract(x, mu),(2*np.pi))==0):
                vm = math.inf
            else:
                vm = 0
    return(vm)

def comp_fn(mu,p_i,x,k,kappa):
    return p_i[k]*vmf_pdf(x,mu[k],kappa[k])
def Ez_fun(mu,p_i,x,k,kappa):
    return comp_fn(mu,p_i,x,k,kappa)/sum(list(map(lambda k:comp_fn(mu,p_i,x,k,kappa) ,list(range(0,knum)) )))
Ez_vals = np.full([n,knum], np.nan) 
for iter in range(1,100):
  # E-step
    for k in range(0,knum):
        Ez_vals[:,k] = list(map(lambda x:Ez_fun(mu,p_i,x,k,kappa) ,dat_ls )) #P(pi of which fn dat came from/parameters)
  
  # M-step
    for k in range(0,knum):
        
        #alpha
        nk = sum(Ez_vals[:,k]) 
        p_i[k] = nk/n
        
        sum_norm = np.linalg.norm(sum(map(lambda x,y:np.dot(x,y),Ez_vals[:,k],dat_vector)))
        mu_vector = sum(map(lambda x,y:np.dot(x,y),Ez_vals[:,k],dat_vector))/sum_norm 
        print("mu_vector",mu_vector)
        mu[k] = math.acos(mu_vector[0])
        print("mean",k,mu[k])
        R = sum_norm / nk
        mle_kappa = (R * (p - R**2)) / (1 - R**2)
        print("kappa",k,mle_kappa)
        kappa[k] = mle_kappa 
        print("aplha",p_i[k])
   
mu[0],mu[1],kappa[0],kappa[1],p_i[0],p_i[1]<jupyter_output><empty_output><jupyter_text>b)  Method 2 - Generalised von mises fischer (Inputs are vectors or any N-dimensional Data)  > Either one of the following 3 cells containing Simulated data is to be run for different kind of input data >> 3 mixtures of 2 features<jupyter_code>c = [[x,y] for x,y in zip(list(np.random.normal(100,20,500)) , list(np.random.normal(50,20,500)))]

d = [[x,y] for x,y in zip(list(np.random.normal(50,10,1000)) , list(np.random.normal(0,10,1000)))]

e = [[x,y] for x,y in zip(list(np.random.normal(50,10,750)) , list(np.random.normal(100,20,750)))]

dat_raw = c+d+e

main_df = pd.DataFrame(dat_raw,columns=['feature-1','feature-2'])<jupyter_output><empty_output><jupyter_text>>> 2 mixtures of 2 features<jupyter_code>c = [[x,y] for x,y in zip(list(np.random.normal(100,20,500)) , list(np.random.normal(50,20,500)))]

d = [[x,y] for x,y in zip(list(np.random.normal(50,10,1000)) , list(np.random.normal(0,10,1000)))]

dat_raw = c+d

main_df = pd.DataFrame(dat_raw,columns=['feature-1','feature-2'])<jupyter_output><empty_output><jupyter_text>>> 2 mixtures of 3 features<jupyter_code>c = [[x,y,z] for x,y,z in zip(list(np.random.normal(100,20,500)) , list(np.random.normal(50,20,500)),list(np.random.normal(50,10,500)))]

d = [[x,y,z] for x,y,z in zip(list(np.random.normal(50,10,1000)) , list(np.random.normal(0,10,1000)),list(np.random.normal(100,20,1000)))]

dat_raw = c+d

main_df = pd.DataFrame(dat_raw,columns=['feature-1','feature-2','feature-3'])
dat_raw
plt.scatter([a[0] for a in dat_raw],[a[1] for a in dat_raw])

def unit_vec_fn(ls):
    """Convert the list of raw n-dimensional data into
    a list of set of unit vectors 
    """    
    magnitude  = pow(sum([a**2 for a in ls]),0.5)
    unit_vector  = [c/magnitude for c in ls]
    return(unit_vector)

dat_ls = [unit_vec_fn(a) for a in dat_raw]

p=len(dat_ls[0])
dat_ls
def _get_vmf_likelihood_term(x, mu, kappa):
    """returns the likelihood back to the function 'vmf_pdf2()' that calculates the pdf"""
    return exp(kappa * np.dot(mu, x))
def _get_vmf_normalization_numerator(p, kappa): 
    """returns a part of the numerator of the pdf back to the function 'vmf_pdf2()' that calculates the pdf"""
    return kappa ** (0.5*p - 1)
def _get_vmf_normalization_denom(p, kappa):
    """returns a part of the denominator of the pdf back to the function 'vmf_pdf2()' that calculates the pdf"""
    return (2 * pi) ** (0.5*p) * iv(0.5*p-1, kappa)
def vmf_pdf2(x, mu, kappa):
    """
    Pdf of the von Mises-Fisher distribution
    Parameters:
        mu: list, location parameter
        kappa: float/int, scale parameter. Large values of kappa corresponds to lower variance
    Returns:
        list, pdf function
    """
    p=len(dat_ls[0])
    likelihood = _get_vmf_likelihood_term(x, mu, kappa)
    normalization_numerator = _get_vmf_normalization_numerator(p, kappa)
    normalization_denominator = _get_vmf_normalization_denom(p, kappa)
    return likelihood * (normalization_numerator / normalization_denominator)

def comp_fn(mu,p_i,x,k,kappa):
    """Numerator of the Expectation step in the algorithm and thiis function
    returns to the call from 'Ez_fun()' """
    return p_i[k]*vmf_pdf2(x,mu[k],kappa[k])

def Ez_fun(mu,p_i,x,k,kappa):
    """Performs the Expectation step in the EM algorithm"""
    return comp_fn(mu,p_i,x,k,kappa)/sum(list(map(lambda k:comp_fn(mu,p_i,x,k,kappa) ,list(range(0,knum)) )))

t1= time.time()
knum = 3

n = len(dat_ls)
mu=[[0.5,0.5],[0.5,0.5],[0.5,0.5]]

kappa=[20,30,30]
p_i = [0.22,0.33,0.45] #[1/2,1/2]

mu_list=[]
pi_list=[]
kappa_list=[]

Ez_vals = np.full([n,knum], np.nan) 
for iter in range(1,500):
    print(kappa)
  # Expectation (E-step)
    for k in range(0,knum):
        """Performs the Expectation step in the EM algorithm for each data set of points"""
        Ez_vals[:,k] = list(map(lambda x:Ez_fun(mu,p_i,x,k,kappa) ,dat_ls )) #P(pi of which fn dat came from/parameters)

  # Maximization Step (M-step)
    for k in range(0,knum):
        
        #alpha
        nk = sum(Ez_vals[:,k]) 
        p_i[k] = nk/n
        # print(Ez_vals[:,k])
#         sum_norm = np.linalg.norm(nk)
        sum_norm = np.linalg.norm(sum(map(lambda x,y:np.dot(x,y),Ez_vals[:,k],dat_ls)))
        mu[k]= sum(map(lambda x,y:np.dot(x,y),Ez_vals[:,k],dat_ls))/sum_norm 

        #print("mean",k,mu[k])
        R = sum_norm / nk
        mle_kappa = (R * (p - R**2)) / (1 - R**2)
        #print("kappa",k,mle_kappa)
        kappa[k] = mle_kappa 
        #print("aplha",p_i[k])

 
    mu_copy = []
    mu_copy = mu.copy()
    mu_list.append(mu_copy)

    pi_copy = []
    pi_copy = p_i.copy()
    pi_list.append(pi_copy)

    kappa_copy = []
    kappa_copy = kappa.copy()
    kappa_list.append(kappa_copy)

    if(iter!=1):
        summ = 0
        for i in range(0,knum):
            summ = summ + sum([abs(a-b) for a,b in zip(np.squeeze(mu_list).tolist()[-1][i],np.squeeze(mu_list).tolist()[-1][i])]) +abs(pi_list[-1][i]-pi_list[-2][i]) +abs(kappa_list[-1][i]-kappa_list[-2][i])
        if(summ < 0.00001):
            itrn = iter
            break   
            
for i in range(0,knum):
    print("mean %s is:"%(i),list(mu[i]))
    print("kappa %s is: %s"%(i,kappa[i]))
    print("proportion %s is: %s"%(i,p_i[i]))  
t2=time.time()
print("time :",t2-t1)<jupyter_output>[20, 30, 30]
[4.59441966280867, 7.284481900183474, 7.284481900183576]
[4.713560651585678, 6.613141841081929, 6.613141841081784]
[4.796138854115853, 7.598493930095134, 7.598493930094797]
[4.898675543141822, 10.757159997927669, 10.75715999792765]
[5.10105242690865, 20.847648521029335, 20.847648521027814]
[5.600289271500151, 54.22380568696903, 54.223805686963274]
[6.782447378678587, 100.51082835010097, 100.51082835013206]
[8.389977337901424, 114.51683902247699, 114.51683902251163]
[9.68993677354458, 109.79274070581023, 109.79274070581023]
[10.541346196594995, 102.35365052155306, 102.35365052156459]
[11.068411496629189, 96.79606742765061, 96.79606742766501]
[11.382047287252604, 93.25180610855757, 93.25180610856522]
[11.56259029670994, 91.12044253267624, 91.12044253267076]
[11.66462899203168, 89.87640601187213, 89.87640601181447]
[11.721877116709216, 89.16265140391775, 89.16265140385138]
[11.753926121153109, 88.75727882661256, 88.75727882657105]
[11.77186208312616, 88.52841736529632, 88.528[...]
{'fraction_non_alphanumeric': 0.17385729370945843, 'fraction_numerical': 0.09752602202443808, 'mean_word_length': 2.758219954648526, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '134281', 'license': 'MIT', 'n_tokens_mistral': 6299, 'n_tokens_neox': 5485, 'n_words': 1090, 'path': 'explanation-&-verification/Von_Mises(&)VMF_mixture.ipynb', 'repo_name': 'maradnusKihtrak/circularAlgos', 'size': '13258'}