Ejercicios Unidad 4
-------------------

Dia 1
-----------------------------

Escribir una funcion piNumber.m que calcule el numero pi por el metodo de
Monte Carlo. En este metodo, se generan aleatoriamente puntos (x,y) en el
cuadrado 0<(x,y)<1 y se cuenta la fraccion de los que cumplen que x^2+y^2<1. 
Esta fraccion debe ser igual al cociente entre la cuarta parte del area 
del circulo de radio unidad (pi/4) y la del cuadrado de lado uno. Interfaz:

function pi = piNumber(n)
% Finds number pi using Monte Carlo method: chooses random points
% in the square 0<(x,y)<1 and counts how many obey x^2+y^2<1. This
% must be the ratio of one fourth the area of the circle of unit
% radius (pi/4) to the area of the square (one).
% Input
%   n  : Number of Monte Carlo trial points
% Output
%   pi : Estimated value of pi
% Usage
%   pi = piNumber(1000)

-----------------------------

Escribir una funcion meanValue.m que devuelva la media de un conjunto de 
valores y pesos, con la interfaz

function meanX = meanValue(x,w)
% Returns the mean value of a set of values x with weights w.
% Input
%   x : vector of values x_i
%   w : optional vector of weights, with same lenght as x
%       If not present, all weights are assumed equal 
% Output
%   meanX = (sum_i x_i * w_i) / (sum_i w_i)
% Usage
%   meanX = meanValue(x,w)

------------------------------

Escribir una funcion standDev.m que devuelva la desviacion estandar 
de un conjunto de valores y pesos, con la interfaz

function deltaX = standDev(x,w)
% Returns the standard deviation of a set of values x and weights w.
% Input
%   x : vector of values x_i
%   w : optional vector of weights, with same lenght as x
%       If not present, all weights are assumed equal 
% Output
%   deltaX = sqrt( (sum_i (x_i - meanX)^2 * w_i) / (sum_i w_i) )
%   where meanX = (sum_i x_i * w_i) / (sum_i w_i)
% Usage
%    deltaX = standDev(x,w)

Nota: el promedio y la desviacion estandar de una serie de n 
valores x_i, con probabilidades p_i=w_i/sum(w_i), se definen como
   meanX = sum(i=1:n) p_i*x_i
   deltaX = sqrt( sum(i=1:n) p_i*(x_i - meanX)^2 )

-----------------------------

Escribir una funcion piNumberAndError.m que calcule el valor y una 
estimacion del error del numero pi. Para estimar el error, se generan 
m valores de pi por el metodo de Monte Carlo y se calcula la desviacion
tipica entre ellos. Interfaz:

function [value,error] = piNumberAndError(n,m)
% Finds number pi using Monte Carlo method: chooses n random points
% in the square 0<(x,y)<1 and counts how many obey x^2+y^2<1. This
% must be the ratio of one fourth the area of the circle of unit
% radius (pi/4) to the area of the square (one). This process is
% repeated m times and the standard deviation of the m results is 
% returned as error.
% Input
%   n  : Number of Monte Carlo trial points in each sub-calculation
%   m  : Number of sub-calculations
% Output
%   value : Estimated value of pi
%   error : Estimated error of pi
% Usage
%   [val,err] = piNumberAndError(1000,10)

-----------------------------

Escribir una funcion coin.m que genere un vector aleatorio de ceros y unos
con la siguiente interfaz

function x = coin(n)
% Returns a vector x of n random bits (zero or one)
% Input
%   n : number of values to be returned
% Output
%   x : row vector of n random bits (zero or one)
% Usage
%   x = coin(100)

Comprobar que, cuando n se hace grande, el promedio y la desviacion estandar
de los vectores devueltos por coin tienden a los mismos valores que 
meanValue(x) y standDev(x), con x=[0,1] (cara=0 y cruz=1 con igual 
probabilidad). 

-----------------------------

Escribir una funcion loadedCoin.m que genere un vector aleatorio de ceros 
y unos con la siguiente interfaz

function x = loadedCoin(w,n)
% Returns a vector x of n random bits (zero or one)
% Input
%   w : weight (probability) of heads result (x=1)
%   n : number of values to be returned
% Output
%   x : row vector of n random bits (zero or one)
% Usage
%   x = loadedCoin(0.55,100)

Comprobar que, cuando n se hace grande, el promedio y la desviacion estandar
de los vectores devueltos por coin tienden a los mismos valores que 
meanValue(x,p) y standDev(x,p), con x=[0,1], p=[1-w,w] (cruz=1 con
probabilidad w).

-----------------------------

Escribir una funcion loadedDice.m que genere un vector aleatorio de 
n tiradas de dados (valores de 1 a 6) con la siguiente interfaz

function x = loadedDice(w,n)
% Returns a vector x of n throws of a loaded dice
% Input
%   w : vector of weights of each result (1 to 6)
%   n : number of values (dice throws) to be returned
% Output
%   x : row vector of n throw results (values 1 to 6)
% Usage
%   x = loadedDice([0.16 0.16 0.16 0.16 0.16 0.20],100)

-------------------------------
Dia 2
-----------------------------

El teorema del limite central o teorema central del limite indica que, 
en condiciones muy generales, si Xn es el promedio de n variables aleatorias 
independientes x, entonces la funcion de distribucion de Xn "se aproxima bien" 
a una distribucion normal (tambien llamada distribucion gaussiana, curva de 
Gauss o campana de Gauss) cuando n es suficientemente grande. La desviacion 
estandar de Xn (la anchura de la curva gaussiana) es igual a la desviacion 
estandar de x dividida la raiz cuadrada de n.

-----------------------------

Escribir una funcion randomAvgs.m que devuelva nAvgs promedios, cada uno de 
nVals valores de la variable x, tomados aleatoriamente con pesos w.

function xAvg = randomAvgs(x,w,nVals,nAvgs)
% Returns a vector of nAvgs averages, each of nVals values, taken randomly
% from a vector x, with relative probabilities w.
% Input
%   x     : vector of possible values
%   w     : vector of weights, with the same lenght as x
%   nVals : number of values in each average
%   nAvgs : number of averages to be calculated (lenght of vector xAvg)
% Output 
%   xAvg  : row vector, of length nAvgs, with random averages
% Usage
%   xAvg = randomAvgs(x,w,nVals,nAvgs)

-----------------------------

Escribir un script plotAvgHistograms.m que, usando la funcion hist de MatLab, 
dibuje histogramas de promedios de bits (x=[0,1],w=[1/2,1/2]) devueltos por 
randomAvgs para nVals = 16, 64, 256 y 1024. Hacer lo mismo para promedios
de tiradas de un dado (x=(1:6),w=1/6).

Comprobar que, para valores suficientemente grandes de nVals y nAvgs, 
la forma de los histogramas es una gausiana
  P(xAvg) = exp(-(xAvg-x0).^2/(2*dxAvg^2)) / sqrt(2*pi*dxAvg^2)
siendo x0 el valor medio de x, y dxAvg la desviacion tipica de la variable 
xAvg, con dxAvg=dx/sqrt(nVals), siendo dx la desviacion tipica de x.

-------------------------------
Dia 3
-----------------------------

La fuerza descendente sobre un cuerpo cayendo por el aire es F = mg - cv. 
El segundo termino se debe al rozamiento con el aire y, para v pequeña, es
mucho menor que el primero. Si lo despreciaramos, la velocidad seria v=~gt.
Si sustituimos esta aproximacion en la expresion de F, la aceleracion sera
a = F/m =~ g - cgt/m. Integrando dos veces obtenemos v =~ gt - (cg/2m)t^2,
x =~ (g/2)t^2 - (cg/6m)t^3.
En una reconstruccion del experimento de Galileo, se deja caer una piedra
de 325 g desde una torre y se mide el tiempo en que pasa por cada piso, 
obteniendose la siguiente tabla

  x (m)    t (s)
 ---------------
   0.0     0.00
   3.2     0.83
   6.4     1.16
   9.6     1.45
  12.8     1.65
  16.0     1.88
 ---------------

Resolver el sistema de ecuaciones x = c2*t^2 + c3*t^3 en funcion de c2 y c3.
Como hay mas ecuaciones que incognitas, usar el operador \ de matlab, que
resuelve el problema como ajuste de minimos cuadrados. Calcular g y c 
a partir de c2 y c3.

Dibujar la curva c2*t^2 + c3*t^3 junto con los puntos (t_i,x_i) de la tabla.
Calcular el error cuadratico medio sqrt(sum((x-(c2*t^2 + c3*t^3))^2) entre
los puntos x de la tabla y su curva de ajuste.

Para estimar el error del ajuste, repetirlo seis veces, excluyendo un punto
cada vez. Calcular el valor medio y la desviacion tipica de los seis valores
de g y c.

En un programa regressFall.m, usar la funcion regress de Matlab para repetir 
el ajuste, obteniendo ademas el intervalo de confianza de c2, c3, g, y c,
asi como las desviaciones tipicas de g y c (la anchura del intervalo de 
confianza del 95% es 3.92 desviaciones tipicas). Dibujar las curvas de 
ajuste para los valores maximos y minimos de c2 y c3.

-----------------------------

Volvemos a considerar ahora la ecuacion original no aproximada
   F = m*g - c*v = m*dv/dt. 
Integrando en t obtenemos
   v = vf + (v0-vf)*exp(-t/tau)
siendo tau=m/c, v0 la velocidad inicial, y vf=g*tau la velocidad en t=inf.
Volviendo a integrar en t obtenemos
   x = x0 + vf*t + (v0-vf)*tau*(1-exp(-t/tau))

a) Comprobar con lapiz y papel que dx/dt=v y dv/dt=a=F/m.

b) Escribir una funcion con la siguiente interfaz
      function x = fall( params, t )
      % Position of a mass falling in a viscuos fluid, with x(0)=v(0)=0
      % Input:
      %   params(2) : parameters g and tau, in m/s^2 and s
      %   t         : time(s) in sec
      % Output:
      %   x         : vertical position(s) in meters (positive downwards)

c) En un programa fitExactFall.m, usar la funcion nlinfit de matlab 
para ajustar los parametros g y c. Dibujar la curva xfit=fall(params,t) 
junto con los puntos (t_i,x_i) de la tabla. Calcular el error cuadratico 
medio sqrt(sum(xfit(t_i)-x_i)^2)) entre los puntos x de la tabla y su curva 
de ajuste.

-----------------------------