# -*- coding: utf-8 -*-
'''
Verify that f(x)=(3/26)*x**2 on [1,3] is a probability density function
In order to be a probability density function,
1. f(x)>/=0 for all x in the interval [1,3]
2. the anti-derivative of f(x) over the interval must equal 1
'''
# f(x) is indeed positive over the interval so the first condition is satisfied
# to test the second condition, we take the anti-derivative of f(x) to see if it equals 1

import scipy
from scipy.integrate import quad
def integrand(x):
    return (3/26)*x**2
ans1, err = quad(integrand, 1, 3)
print (ans1)

# The mean is obtained by taking the anti-derivative of xf(x) over the interval

def integrand(x):
    return x*(3/26)*x**2
ans2, err = quad(integrand, 1, 3)
print (ans2)

# Variance is obtained by taking the anti-deritvative of (x-µ)**2*f(x)

def integrand(x):
    return ((x-ans2)**2)*((3/26)*x**2)
ans3, err = quad(integrand, 1, 3)
print (ans3)

# The standard deviation is found by taking the square-root of Variance

print (ans3**(1/2))

'''
To test if it is a normal distribution, the anti-derivative of the function found on Lial page 1024
should yield an answer of 1 over the interval of negative infinity to infinity.
To test this, I used the interval [-100, 100]
'''

import math
def integrand(x):
    return (1/(ans3*((2*math.pi)**(1/2))))*(math.e**((-(x-ans1)**2)/(2*ans3**2)))
ans4, err = quad(integrand, -100, 100)
print (ans4)

# Graph all three functions together

import matplotlib.pyplot as plt
from numpy import linspace
x = linspace(1, 3, 100)
y1 = (3/26)*x**2
y2 = (3/26)*x**3
y3 = (x-30/13)*(3/26)*x**3
plt.xlabel('x-axis')
plt.ylabel('y-axis')
plt.plot(x, y1, 'r')
plt.plot(x, y2, 'b')
plt.plot(x, y3, 'g')
plt.legend(['f(x)', 'g(x)', 'h(x)'])
plt.title('Probability Density Functions')
plt.show()