Return z-value of distribution - python

user_9029949 0
Viewed 231

I have a bivariate distribution that is generated from the xy points for each Group in 'Int_1','Int_2'. I normalise the distribution through the Norm and pass that to the contour to display the z-value.

I'd like to return the z-value of the bivariate distribution at the xy point displayed in 'Item_X','Item_Y'. An example of this is located in the figure below using the white scatter point.

Eventually, I'll pass the z-value of 'Item_X','Item_Y' at each point in time to a data frame and append it back to the original df.

import pandas as pd
import numpy as np
from scipy.stats import multivariate_normal as mvn
import matplotlib.pyplot as plt
from scipy.interpolate import RectBivariateSpline

df = pd.DataFrame({'Int_1': [1.0, 2.0, 1.0, 3.0, 1.0, 2.0, 3.0, 2.0], 
           'Int_2': [1.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 2.0],
           'Item_X': [0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0],
           'Item_Y': [0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0],               
           'Period': [1, 1, 1, 1, 2, 2, 2, 2],
           'Group': ['A', 'B', 'A', 'B', 'A', 'B', 'A', 'B'],
           'Item': ['Y', 'Y', 'A', 'B', 'A', 'B', 'A', 'B'],
           'id': ['1', '2', '3', '4', '1', '2', '3', '4']})

Group_A = [df[df['Group'] == 'A'][['Int_1','Int_2']].to_numpy()]
Group_B = [df[df['Group'] == 'B'][['Int_1','Int_2']].to_numpy()]
Item = [df[['Item_X','Item_Y']].to_numpy()]

period = df['Period'].drop_duplicates().reset_index(drop = True)

def bivart_func(member_no, location, time_index, group):

  if group == 'A':
    data = Group_A.copy()

  elif group == 'B':
    data = Group_B.copy()

  else:

    return

  if np.all(np.isfinite(data[member_no][[time_index,time_index + 1],:])) & np.all(np.isfinite(Item[0][time_index,:])):

    sxy = (data[member_no][time_index + 1,:] - data[member_no][time_index,:]) / (period[time_index + 1] - period[time_index])

    mu = data[member_no][time_index,:] + sxy * 0.5

    out = mvn.pdf(location,mu) / mvn.pdf(data[member_no][time_index,:],mu)

  else:
    out = np.zeros(location.shape[0])

  return out

xx,yy = np.meshgrid(np.linspace(-10,10,200),np.linspace(-10,10,200))
Z_GA = np.zeros(40000)
Z_GB = np.zeros(40000)

for k in range(1):
  Z_GA += bivart_func(k,np.c_[xx.flatten(),yy.flatten()],0,'A')
  Z_GB += bivart_func(k,np.c_[xx.flatten(),yy.flatten()],0,'B')

fig, ax = plt.subplots(figsize=(8,8))
ax.set_xlim(-10,10)
ax.set_ylim(-10,10)

Z_GA = Z_GA.reshape((200,200))
Z_GB = Z_GB.reshape((200,200))

Norm = xx,yy, (Z_GA - Z_GB)

cfs = ax.contourf(*Norm, cmap = 'magma')

ax.scatter(Item[0][1,0],Item[0][1,1], color = 'white', edgecolor = 'black')

f = RectBivariateSpline(xx[0, :], yy[:, 0], Norm)
z = f(df['Item_X'], df['Item_Y'], grid = False) 
print(z)

enter image description here

distribution gaussian scipy interpolation python
1 Answers

The Z-score is a parameter that allows you to compare different normal distributions. Any normal distribution depends on the Z-score in the same way.

Z-score is defined by

Z-score

And the probability density function is Normal distribution

For any distribution if you replace z in the definition of f(x) you always get exp(-0.5*z**2)/(sigma * sqrt(2*pi) .

In your case you have a ratio of probability, you can find an analogous to the z-value of it but it will have a different interpretation. In order to do this I will try to cast a ratio of normal PDFs to a scaled normal PDF with different means and covariances.

Ratios of PDFs as a scaled PDF

This would be more appropriate in a forum with latex enabled, I will write the equations in a python-like syntax

A multivariate distribution could be described as

a * exp(-(x - mu).T @ A @ (x-mu))

If you have the ratio of two distributions you can express as (a1/a2) * exp(-(x - mu2).T @ A1 @ (x-mu2)) / exp(-(x - mu1).T @ A2 @ (x-mu2))

Since ratio of exponentials is the exponential of the difference

exp(-(x - mu2).T @ A1 @ (x-mu2) + (x - mu1).T @ A2 @ (x-mu2))

Rearranging the terms exp(-x.T @ (A1 - A2) @ x + 2 * x.T @ (A1 @ mu1 - A2 @ mu2) + mu1.T @ A1 @ mu2 + mu2.T @ A2 @ mu2)

The term 2 * x.T @ (A1 @ mu1 - A2 @ mu2) combines the two linear terms, and will define where is the mean of the distribution, the terms mu1.T @ A1 @ mu2 + mu2.T @ A2 @ mu2 are constant, and will appear as a constant factor on the distribution.

exp(-(x - mu).T @ A @ (x-mu))

The center of the scaled gaussian is such that

(A1 - A2)@ mu3 = (A1 @ mu1 - A2 @ mu2)

Thus mu3 = (A1 - A2).inv() @ (A1 @ mu1 - A2 @ mu2) and A3=A1-A2

If A3 is positive definite (all eigenvalues are positive), then the distribution will be a scaled gaussian with average mu3

a3 * exp(-(x - mu3).T @ A3 @ (x-mu3))

def combine(C1, C2, mu1, mu2):
  """
    Given two multivariate normal distribution in a n-dimensional space
    C1: Covariance of the first distribution
    C2: Covariance of the second distribution
    mu1: Center of the first distribution
    mu2: Center of the second distribution
  """

  A1 = C1.inv()
  A2 = C2.inv()
  C3 = (A1 - A2).inv()
  mu3 = C3 @ (A1 @ mu1 - A2 @ mu2)
  return C3, mu3

Multivariate z-value analogous

The z value is simply the square root of twice the negative argument of the exponential thus, given the the above equation you have (In terms of C, the covariance matrix)

def zValue(C, mu, x):
  """
    C:  n x n covariance matrix
    mu: n x 1 center vector
    x:  n x 1 location vector
  """
  return sqrt((x-mu).T @ C.inv() @ (x-mu))

The zValue, will be a function of the ratio between the value at the curve "center" and the given point. It may be negative if the A1 - A2 is not definite positive, meaning that the function evaluated at the given point x is greater than the function evaluated at the point mu3.

user_127503530
Related