#!/usr/bin/env python
# encoding: utf-8
"""
untitled.py

Created by Austin Holland on 2012-07-31.
Copyright (c) 2012 __MyCompanyName__. All rights reserved.
"""

import sys
import os
import numpy as np

def weighted_moment(m,m_sigma,b=1.0):
  N=len(m)
  beta=b*np.log(10)
  m_varscaled=(m_sigma**2)-np.var(m)
  mx_var=np.sum(m_varscaled)
  em_sum=0.
  em_var=0.
  for i in np.arange(0,N):
    em_var+=1./(m_varscaled[i])  

  em_var=1./(em_var)

  for i in np.arange(0,N):
    em_sum+=(em_var/m_varscaled[i])*m[i]
 
  em=em_sum-(N-1)*beta*em_var
  
  return em,np.sqrt(em_var)
  
def myweighted_moment(m,m_sigma,b=1.0):
  N=len(m)
  beta=b*np.log(10)
  numsum=0.
  denomsum=0.
  m_var=(m_sigma**2)
  for i in np.arange(0,N):
    numsum+=m[i]/m_var[i]
    denomsum+=1/m_var[i]
    
  me=numsum/denomsum
  m_uncert=np.sqrt(1./denomsum)
  cor=(N-1)*beta*np.var(m)
  return me,m_uncert,cor
  
def main():
  pass


if __name__ == '__main__':
  m=np.array([4.8,5.0,4.6])
  m_sigma=np.array([0.25,0.25,0.5])
  print m,m_sigma
  print weighted_moment(m,m_sigma)
  print "My moment"
  print myweighted_moment(m,m_sigma)
  
  m=np.array([1.584,1.684,2.3168]) 
  m_sigma=np.array([0.24,0.24,0.25])
  print "\n"
  print m,m_sigma
#  print weighted_moment(m,m_sigma)
  print myweighted_moment(m,m_sigma)