lab2.py

import numpy as np
import math

n = 400


def initX(n, mu, sigma):
    return np.random.multivariate_normal(mu, sigma, n)


def initR(n, size):
    res = []
    for i in xrange(0, n):
        row = []
        for j in xrange(0, size):
            val = np.random.random_sample()
            if val <= 0.5:
                row.append(1.0)  # obs
            else:
                row.append(0.0)
        res.append(row)

    return res


def complexMethod(x, r):
    x_obs = []

    for i, r_i in enumerate(r):
        if sum(r_i) == len(r_i):
            x_obs.append(x[i])

    x_obs_mean = []
    sums = sum(x_obs[:])
    for i in xrange(len(x[0])):
        x_obs_mean.append((sums[i]) / len(x_obs))

    print "Mean = " + str(x_obs_mean)

    x_obs_cov = []
    for i in xrange(len(x_obs[0])):
        row = []
        for j in xrange(len(x_obs[0])):
            val = 0
            for k in xrange(len(x_obs)):
                val += (x_obs[k][i] - x_obs_mean[i]) * (x_obs[k][j] - x_obs_mean[j])
            row.append(val / float(len(x) - 1))

        x_obs_cov.append(row)

    print "Covariance = " + str(x_obs_cov)

    x_obs_corr = []
    for i in xrange(len(x[0])):
        row = []
        for j in xrange(len(x[0])):
            row.append(x_obs_cov[i][j] / math.sqrt(x_obs_cov[i][i] * x_obs_cov[j][j]))
        x_obs_corr.append(row)

    print "Correlation = " + str(x_obs_corr)
    return x_obs_mean, x_obs_cov, x_obs_corr


def availableMethod(x, r, printRes=True):
    x_mean = []
    for i in xrange(len(x[0])):
        numerator = 0
        denominator = 0
        for k in xrange(len(x)):
            numerator += r[k][i] * x[k][i]
            denominator += r[k][i]
        x_mean.append(numerator / float(denominator))

    if printRes:
        print "Mean = " + str(x_mean)

    x_cov = []
    for i in xrange(len(x[0])):
        row = []
        for j in xrange(len(x[0])):
            val = 0
            den = 0
            for k in xrange(len(x)):
                val += (r[k][i] * r[k][j]) * (x[k][i] - x_mean[i]) * (x[k][j] - x_mean[j])
                den += (r[k][i] * r[k][j])
            row.append(val / float(den - 1))

        x_cov.append(row)

    if printRes:
        print "Covariance = " + str(x_cov)

    x_corr = []
    for i in xrange(len(x[0])):
        row = []
        for j in xrange(len(x[0])):
            row.append(x_cov[i][j] / math.sqrt(x_cov[i][i] * x_cov[j][j]))
        x_corr.append(row)

    if printRes:
        print "Correlation = " + str(x_corr)

    return x_mean, x_cov, x_corr


def f1(x, r, x_mean, i, j, k, l):
    res = 0
    for n in xrange(len(x)):
        res += (r[n][k] * r[n][l] * (x[n][i] - x_mean[i][i]) * (x[n][j] - x_mean[j][j]))
    return res


def f2(r, k, l):
    res = 0
    for i in xrange(len(r)):
        res += r[i][k] * r[i][l]
    return res - 1


# todort
def availableMethod2(x, r):
    x_mean = []
    for i in xrange(len(x[0])):
        row = []
        for j in xrange(len(x[0])):
            numerator = 0
            denominator = 0
            for k in xrange(len(x)):
                numerator += r[k][i] * r[k][j] * x[k][i]
                denominator += r[k][i] * r[k][j]
            row.append(numerator / float(denominator))
        x_mean.append(row)

    print "Mean = " + str(x_mean)

    x_cov = []
    for i in xrange(len(x[0])):
        row = []
        for j in xrange(len(x[0])):
            val = 0
            den = 0
            for k in xrange(len(x)):
                val += (r[k][i] * r[k][j]) * (x[k][i] - x_mean[i][j]) * (x[k][j] - x_mean[j][i])
                den += (r[k][i] * r[k][j])
            row.append(val / float(den - 1))

        x_cov.append(row)

    print "Covariance wave = " + str(x_cov)

    x_corr_wave = []
    for i in xrange(len(x[0])):
        row = []
        for j in xrange(len(x[0])):
            row.append(x_cov[i][j] / math.sqrt(x_cov[i][i] * x_cov[j][j]))
        x_corr_wave.append(row)

    print "Correlation wave = " + str(x_corr_wave)

    x_corr_star_wave = []
    for i in xrange(len(x[0])):
        row = []
        for j in xrange(len(x[0])):
            row.append(x_cov[i][j] / math.sqrt(
                (f1(x, r, x_mean, i, i, i, j) / float(f2(r, i, j))) * (
                    f1(x, r, x_mean, j, j, i, j) / float(f2(r, i, j)))))
        x_corr_star_wave.append(row)

    print "Correlation star wave = " + str(x_corr_star_wave)

    x_corr_star = []
    x_mean_available, x_cov_available, x_corr_available = availableMethod(x, r, False)
    for i in xrange(len(x[0])):
        row = []
        for j in xrange(len(x[0])):
            row.append(x_cov_available[i][j] / math.sqrt(
                (f1(x, r, x_mean, i, i, i, j) / float(f2(r, i, j))) * (
                    f1(x, r, x_mean, j, j, i, j) / float(f2(r, i, j)))))
        x_corr_star.append(row)

    print "Correlation star = " + str(x_corr_star)


def monteCarlo(n, mean, var):
    list_mean = []
    list_var = []
    for i in xrange(n):
        x = initX(n, mean, var)
        r = initR(n, 2)
        x_mean, x_var, corr = availableMethod(x, r)
        list_mean.append(x_mean)
        list_var.append(x_var)

    return sum([np.array(np.array(list_mean[i]) - np.array(mean)) for i in xrange(n)]) / float(n), sum(
        [np.array(np.array(list_var[i]) - np.array(mean)) for i in xrange(n)]) / float(n)


def main():
    x = initX(n, [1.0, -1.0], [[0.66, 0.1], [0.1, 0.33]])
    r = initR(n, 2)
    print "complex"
    complexMethod(x, r)
    print "available"
    availableMethod(x, r)
    print "complex"
    availableMethod2(x, r)

    # x = initX(n, [1.0, 1.0, 1.0], [
    #     [1.0, 0.1, 0.2],
    #     [0.1, 1.0, 0.3],
    #     [0.2, 0.3, 1.0]
    # ])
    # r = initR(n, 3)
    # print "complex"
    # complexMethod(x, r)
    # print "available"
    # availableMethod(x, r)
    # print "complex"
    # availableMethod2(x, r)

    # print "monte-carlo:"
    # _m, _v = monteCarlo(n, [1.0, 1.0], [[1.0, 1.0], [1.0, 1.0]])
    # print "mean = " + str(_m)
    # print "covar = " + str(_v)


main()