# -*- coding: utf-8 -*-
"""James-Stein-Estimator.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1UjZf01Tnf3DW8ZjpP6Osl0LDk_jLCFTP
"""

import numpy as np
import random
import matplotlib.pyplot as plt

def pos(x):
  if x>0:
    return x
  else:
    return 0

def risk_comparison(k, n):

    # Intiliazing vector of means centered at 0
    mu = [0 for i in range(k)]

    # generate 1000 samples
    X = np.random.multivariate_normal(mu, np.identity(k), n)

    # MLE estimator and Risk
    theta_mle = np.mean(X, axis=0)
    risk_mle = np.sum(np.square(mu - theta_mle))
    # print ("Risk MLE: ", risk_mle)

    # James Stein estimator and risk
    theta_mle_norm = np.linalg.norm(theta_mle, ord=2)
    theta_js = pos((1 - ((k-2)/(n*(np.square(theta_mle_norm)))))) * theta_mle
    risk_js = np.sum(np.square(mu - theta_js))
    # print ("Risk James-Stein", risk_js)

    # print ("Diff in risk: ", risk_js - risk_mle)
    return risk_js, risk_mle

# Testing
k = 5
n = 10
risk_comparison(k, n)

# Function to plot for k=1 to 100

def risk():
  n = 1000        # Number of samples for each k
  js_risk = []
  mle_risk = []
  for k in range(1,100):
    js, mle = risk_comparison(k, n)
    js_risk.append(js)
    mle_risk.append(mle)

  # Plotting the risk estimates
  plt.figure(figsize=(8,8))
  plt.plot(js_risk,label='James-stein risk')
  plt.plot(mle_risk,label='MLE risk')
  plt.legend(('James-stein risk', 'MLE risk'),
           loc='upper center')
  plt.xlabel('K parameter theta')
  plt.title('Comparison of James-stein estimator and MLE estimator')
  plt.ylabel('Risk values')
  plt.show()

risk()


# -*- coding: utf-8 -*-
"""James-Stein-Estimator.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1UjZf01Tnf3DW8ZjpP6Osl0LDk_jLCFTP
"""

import numpy as np
import random
import matplotlib.pyplot as plt

def pos(x):
  if x>0:
    return x
  else:
    return 0

def risk_comparison(k, n):

    # Intiliazing vector of means centered at 0
    mu = [0 for i in range(k)]

    # generate 1000 samples
    X = np.random.multivariate_normal(mu, np.identity(k), n)

    # MLE estimator and Risk
    theta_mle = np.mean(X, axis=0)
    risk_mle = np.sum(np.square(mu - theta_mle))
    # print ("Risk MLE: ", risk_mle)

    # James Stein estimator and risk
    theta_mle_norm = np.linalg.norm(theta_mle, ord=2)
    theta_js = pos((1 - ((k-2)/(n*(np.square(theta_mle_norm)))))) * theta_mle
    risk_js = np.sum(np.square(mu - theta_js))
    # print ("Risk James-Stein", risk_js)

    # print ("Diff in risk: ", risk_js - risk_mle)
    return risk_js, risk_mle

# Testing
k = 5
n = 10
risk_comparison(k, n)

# Function to plot for k=1 to 100

def risk():
  n = 1000        # Number of samples for each k
  js_risk = []
  mle_risk = []
  for k in range(1,100):
    js, mle = risk_comparison(k, n)
    js_risk.append(js)
    mle_risk.append(mle)

  # Plotting the risk estimates
  plt.figure(figsize=(8,8))
  plt.plot(js_risk,label='James-stein risk')
  plt.plot(mle_risk,label='MLE risk')
  plt.legend(('James-stein risk', 'MLE risk'),
           loc='upper center')
  plt.xlabel('K parameter theta')
  plt.title('Comparison of James-stein estimator and MLE estimator')
  plt.ylabel('Risk values')
  plt.show()

risk()


# -*- coding: utf-8 -*- James-Stein-Estimator.ipynb  Automatically generated by Colaboratory.  Original file is located at     https://colab.research.google.com/drive/1UjZf01Tnf3DW8ZjpP6Osl0LDk_jLCFTP # Intiliazing vector of means centered at 0 # generate 1000 samples # MLE estimator and Risk # print ("Risk MLE: ", risk_mle) # James Stein estimator and risk # print ("Risk James-Stein", risk_js) # print ("Diff in risk: ", risk_js - risk_mle) # Testing # Function to plot for k=1 to 100 # Number of samples for each k # Plotting the risk estimates