Microeconomics 05

From the previous discussion, we can easily see that if unconstrained, consumers will simply buy an infinite amount of goods (by choosing the highest indifference curve).

However, due to a limited budget, we will have to choose the best combination that will maximize our budget and utility. This is called Constraint Maximization.

Budget Constraints

Typically, we can represent  a constraint as:

Where P is the number of pizzas, and Pp is the price of each pizza.

Suppose that a pizza costs $20, and a movie costs $10, and we have $100 to spend, then we will have a graph as such:

We will be drawing $100 = $20P + $10M. The slope of this line is -Pm/Pp.

The price ratio is -0.5, which is the Marginal Rate of Transformation (MRT). It is the rate at which we can transform pizzas into movies.

The Opportunity Cost of something is the value of foregone alternative. For example, if we forgo a pizza, then we are forgoing 2 movies. Or from the other perspective, a movie costs half a pizza. This is the Opportunity Cost of the situation. Opportunity costs only make sense if we have a limited budget.

Suppose that the price of pizza rose to $30. We will have the following graph now:

With the $100, we have things to choose from now. In other words, our Opportunity Set has become more restricted.

We can also have a restricted Opportunity Set if our budget drops. Keeping the price constant as the first example, we will have the following graph if we drop our budget to $80. In this case, the slope remains the same.

We want to find out what is the furthest indifference curve given a budget constraint. An example when we combine the indifference curves and the budget constraint looks like this.

(Note that the x-axis represents movies, while y-axis represents pizza).

We can afford anything under our budget line. Notice that there are three points on our budget line itself. Points A and C are equivalent, but because we want the farthest indifference line (due to Non-Satiation), point D is actually the best we can get.

The best indifference curve we can get is one where the budget constraint line is tangent to the curve. In order to perform utility maximization, we need to get these things to be tangent by ensuring that the MRS is equal to the MRT. By setting MRS equal to MRT, we are setting marginal benefits equal to marginal costs.

Deriving Marginal Utilities, MRS and MRT

The marginal at each point is the partial derivative of the functional utility with respect to the item we are interested in. For example, if we are interested in the marginal utility function of pizza then we differentiate sqrt(PM) with respect to P.

Points A and C are both undesirable. It is clear that the MRT for the entire budget line is -0.5. Let's work out the MRS for all the points:

An MRS of -2 means that we are willing to trade 2 pizzas for a movie. However, the market MRT of -0.5 only requires us to trade half a pizza for a movie.

There are times where we will end up with negative quantities when we solve for the optimum point. It may be that we are looking at a problem with a corner solution.

Remember that MRS is the rate the consumer is willing to trade the y-axis for the x-axis. On the other hand, MRT is the actual market value ratio between the y-axis and the x-axis.

Note that B, D and E both satisfies that MRT equals MRS, but we want the farthest indifference point that is still within our budget. That is the basic goal of utility maximization.

Microeconomics 04

Consumer Theory

In the previous lectures, we simply made use of given supply and demand curves. Now we are going to start looking at how the curves are made. We shall start by focusing on the demand curves in the study of consumer theory.

All consumer-related behavior in economics revolves around utility maximization. In order to do these, we need to take into account consumer preferences, budget constraints, and perform constraint optimization. In the most basic sense, we will look at the tradeoff between two goods and how the consumer will choose between them.

In summary, we will need to:
1) Make assumptions about preferences
2) Translate these assumptions into utility functions
3) Add budget constraints and perform maximization

In order to model preferences across goods, we will need to impose the following assumptions to make our modeling simpler:
1) Completeness - When we compare two bundles of goods, we would prefer one or the other but never equally. In other words, you can always make a choice between them.
2) Transitivity - If we prefer x to y, and y to z, then we will prefer x to z.
3) Non-Satiation - More is always better and we will never be satisfied. We will never turn down having more.

Indifference Curves

These curves are somewhat like preference maps. These are the graphical representation of people's preferences.

Suppose that your parents gave you some money and you want to decide between pizza or movies. We can have the following choices.

Assume that we are indifferent between two pizzas and one movie, and one pizza and two movies. Clearly we would prefer two pizzas and two movies better than either of them. An indifference curve is a curve which shows all combinations of assumptions that the consumer is indifferent. Example curves based on the above statements can be shown as such:

From the assumptions of consumers above, we have the following properties of indifference curves:
1) Due to the non-satiation assumption, consumers will prefer higher indifference curves.
2) Due to non-satiation, we can also say that indifference curves are always downward sloping. If it is upward sloping, then there can be a case where we are indifferent between 2 pizza and 2 movies and 1 pizza and 2 movies.
3) Indifference curves cannot cross. If the curves cross, there would be a case where we are indifferent between two choices even though one of it has more, which violates non-satiation.
4) Completeness also implies that we cannot have more than one indifference curves through a point.

Utility

Utility functions are just a mathematical representation of the preference maps. We simply need to maximize the utility function to see what the user would choose.

An example utility function can be U = sqrt(M*P). This is something that is empirical that we come up with and this is not the only utility function we can come up with. However, this is consistent with the indifference maps.

Marginal utility refers to the amount utility changes with each change of the unit. In other words, it is the (partial) derivative of the utility function.

An example of a diminishing marginal utility is shown below:

We can see that each additional movie improves the utility, but it increases at a diminishing rate. Marginal utility would usually be diminishing, but at different rates for different goods. However, marginal utility will always be positive due to non-satiation.

The Delta Marginal Utility graph can show us the actual contributions of each pizza.

Marginal Rate of Substitution

The MRS links the utility to the preference map. The MRS is the slope of the indifference curve, which is dP/dM. It is the rate at which we are willing to trade off the y-axis for the x-axis. For example, it is how many pizzas we are willing to trade off to get another movie. It purely comes from the preferences.

We will refer to the following diagram and attempt to compute the MRS for each segment.

From the first segment to the second, we have an MRS of -2, because we are willing to give up 2 movies for 1 pizzas. However, the MRS of the second to the third segment is -0.5. This is because when we have 4 pizzas, the last pizza has a very low marginal utility. However its, marginal utility is increasing the lower number of pizzas we have.

In general, we have:

Marginal utility is a negative function of quantity. The lower quantity we have, the higher the marginal utility, which is why we flip the X and Y axis.

Microeconomics 03

Elasticity

The magnitude of change we can expect from a supply or demand shock is determined by shape of the two curves.

The elasticity of the curves can determine how much the quantity or price changes upon a shock.

A perfectly inelastic demand curve is shown below.

When there are no substitutes for a particular product, for example body parts for a transplant, then the demand curve for that product will be inelastic. When there is a supply shock, the quantity doesn't change - only the price changes.

On the other hand, we have a perfectly elastic demand curve below.

A very substitutable product will be elastic. A possible example of an elastic good is toilet paper, where we'll simply switch to another brand if a particular brand becomes expensive. Once the price increases, the quantity sold will drop drastically.

Elasticity can be expressed as:

If the quantity falls by 2% for every 1% increase in price, then we have E = -2. Therefore, an inelastic good will have E = 0, and a perfectly elastic good will have E = -infinity. E will typically be between 0 and -infinity.

Suppose now that a producer sells Q goods at price P, they will make revenue R = PQ. The change in revenue with respect to price, would be:

From this formula, if we are a producer, we should only raise prices when E is between 0 and -1.

Estimating Elasticities

Theoretical economics can show us the direction of change, but empirical economics is what we use to show us the absolute values of change. However, the most difficult part of empirical economics is causation vs correlation.

Take the example in the previous article, we saw that the price of pork rose when the price of beef rose.

If we wanted to calculate E here, we might try to plug in the numbers in the formulae. Since both the price and quantity increased in this case, we may say that we have a positive elasticity - higher price led to higher quantity. However, this is one of the most dangerous mistakes we can make. In this case, the increase in demand is actually what's causing the price to increase, and not the other round. We can only use this formula to find the amount of change in quantity per percentage change in price. Since quantity is driving the price in this case, and not the other round, we cannot apply the elasticity formula.

What we are measuring here is actually the elasticity of the supply.

However, when we change the supply as in the following diagram, we can actually get the correct answer.

What we want is to be able to measure the slope of the demand curve. By shifting the supply curve, keeping the demand curve unchanged, we can measure the elasticity.

Taxation Example

Suppose that we have vendors selling 100 million units of a particular type of product for $10 each. The government comes along and imposes a tax on the vendors, asking for $1 per unit. This increase in price effectively shifts the supply curve up. In order to make the same amount, the producers would have to charge $11 for each unit. However, this causes the equilibrium to drop to 97 million units, at $10.5 per unit.

In this case, we can clearly see that the equilibrium points have traveled on the demand curve, so we can calculate the elasticity using the formula.

Since we see that the elasticity is -0.6, this means that we can continue raising prices further to increase profits.

Of course, we are assuming elasticity is constant. However, since elasticity is a curve for straight demand lines, we are actually just measuring the LOCAL elasticity around that price change.

It can be seen that the amount of money the government earns is the final equilibrium quantity multiplied by the amount of tax. This can be represented by the shaded region. We can clearly see that the amount of money the government makes depends on the elasticity of demand.

If demand is perfectly inelastic, then the amount the government makes is simply Q * tax per unit. If the product demand is elastic, then the government will make a lot less money as the price needs to be fixed.

Microeconomics 02

Introduction to Supply and Demand

The twin engines of economics are Supply and Demand. Demand refers to how much people wants something, and Supply refers to how much something is available. Adam Smith first came up with the Supply and Demand model, and when describing his model, he came up with the Water and Diamond paradox.

Water and Diamond Paradox

We can clearly see that water is much more essential to life than diamond, yet diamonds are worth so much more than water. The problem here is that even though there is great demand for water, the supply is much higher than the demand, hence its low price. On the other hand, even though there is much less demand for diamond than water, its supply is much lower than the demand, hence the high price.

Supply and Demand Equilibrium

We can visualize supply and demand on a graph. An example of a supply and demand curve is shown below.

Generic Supply and Demand Graph

The blue line represents the demand curve. The demand curve represents the consumers' willingness to pay for a certain good. As the price increases, the quantity the consumer obtains is likely to decrease.

On the other hand, the red line represents the supply curve. The supply curve represents the willingness of a producer to supply a good. As the price rises, the producers would be more willing to produce more.

When the two lines meet, we have a Supply and Demand Equilibrium. This is where both the suppliers and consumers will be happy.

Equilibrium Shift

Suppose that the graph represents the supply of chicken. Suppose that the supply of pork goes down due to a disease, leading to an increase in pork prices. Since pork and chicken are substitutes for each other, we can expect the demand for chicken to increase.

Increase in Demand

As the demand increases, the demand curve shifts upwards (or outwards). With the supply being fixed, this would mean that suppliers can start charging more as consumers are now more willing to pay for it.

On the other hand, if the supply for pork decreases with the demand being fixed, we can also end up with a higher price as shown in the following graph.

Decrease in Supply

Although prices increased in both cases, the quantity sold differs between them. In order to determine if it's a supply or demand shift, we need to be told the price, as well as the quantity.

Constraints

In certain countries, there is the concept of minimum wage. We can analyze employment using the Supply and Demand model as well.

Constraint with Excess Supply

In this case, the suppliers would the citizens (as they provide man-hours), while the consumers are the firms. In a country without minimum wage, we would expect the equilibrium to follow the supply and demand curves.

However, suppose that we add the constraint of minimum wage. When there is minimum wage, we would be in a state of disequilibrium. Due to the minimum wage, workers would be more willing to work. However, due to the higher costs, firms will be less willing to hire. This would lead to an excess in supply i.e. unemployment.

Notice that the new equilibrium e2 is actually on the demand curve instead of the supply curve. This is because even though there are more workers willing to work, the firms are the ones deciding how many they want to hire.

If, however, suppose that the minimum wage is below the equilibrium wage. In this case, the wages will actually stabilize at the equilibrium e1 instead of at the minimum wage.

In the above example, we talked about excess supply. There can also be a case where we have excess demand. Suppose that we try to model the Supply and Demand curves for oil.

Excess Demand

Initially, we are at the equilibrium point e1. Suppose that there is a worldwide shortage of oil. If the supply for oil decreases, prices are expected to increase. Shifting the supply curve upwards (to indicate the decrease of supply), we will end up at point e2. However, if the government decides to put a cap on the maximum price for gas (limiting it to the original price), the suppliers will now be unwilling to supply as much as before. We now end up with the equilibrium on the supply curve at e3.

Market Efficiency

In general, equilibrium points always results in the highest efficiency in the market. Constraints, therefore, lead to inefficiencies. The efficiency loss in the wage scenario is that even though workers are willing to work at a lower wage, they are now unemployed because of the minimum wage. Another example is where suppliers would be willing to supply more gas (and consumers will be willing to purchase them) but the gas price cap results in the shortage of gas.

A perfectly competitive market refers to a market where producers offer their goods to a wide range of consumers who bid up the price until the highest bidder gets it. An example of a perfectly competitive market is eBay auctions.

Without constraints, i.e. in a perfectly competitive market, the mechanism used for allocation is price. When prices are allowed to swing unconstrained, we will end up at the natural equilibrium. Looking at the gas example, price is no longer the allocation mechanism due to the constraint which leads to gas shortage. People will have to queue up for gas no matter how much money they are willing to pay. When price is not used as the allocation mechanism, there will be more inefficient mechanisms like the queue mechanism.

There is always a trade-off between efficiency and equity. Even though the market is most efficient at the equilibrium, but it may not be completely fair. For example, without a minimum wage, there are people who could be exploited. When equity comes into play, things become very complicated. For example, in order to shift the supply curve downwards, the government can provide subsidies to the oil companies. However, this would mean that people will be taxed more heavily on the other end.

Summary

In summary, the market is always going to attempt to reach equilibrium whenever they can. Whenever a constraint causes a disequilibrium, the component that is lacking (either supply or demand) will determine where the new equilibrium will be. Equilibrium points are points of highest market efficiency, and constraints, in general, lead to inefficiencies. However, inefficiency may not necessarily be bad as it leads to equity.

Microeconomics 01

The study of microeconomics involves constrained optimization in face of scarcity. We build models on how consumers and producers behave in order to study their relationships. However, unlike engineering models, these models are never precise, and can only capture the main insights and never the minute details.

The main constraint faced by consumers is their limited budget, which we will optimize through utility maximization. We want to maximize their utility subject to a budget constraint. Firms on the other hand focus on maximizing profits subject to both consumer demands and input costs.

In Microeconomics, prices play three fundamental roles:
1) What goods and services should be produced?
2) How do we produce these goods and services?
3) Who will make use of these goods and services?

The above three question can be summarized into the term "allocation". As consumers and firms interact in the marketplace, these questions will lead to a set of prices which best satisfies both the consumers and producers.

If we want to create a product, we need to analyze the consumer's reception towards it. Next, we will need to look at all the input costs to see if we will be able to set a price that allows us to profit, while making it affordable to the target market.

When dealing with economics, we have to face both the theoretical and empirical sides. Theoretical economics involves building models to explain the world, while empirical economics involves testing those models to judge its accuracy. That would mean that, of course, models built in theoretical economics should be testable.

Economics can also be positive or normative. Positive describes what the world are, and normative describes what the world should be.

We can almost model any kind of decision. An example model for a decision-making process for buying a new vs used product can be:
1) What are your personal preferences?
2) Are you risk-averse or risk-loving?
3) Is your budget limited or unlimited?

By comparing the price of the new and used product, we weigh that difference in price with the above questions to result in an optimization problem.

Even though we do not really go through these thought processes, but we behave AS IF we have solved these optimization problems. This is known as Milton Friedman's "As If" principle.

Physics 01

Introduction to Dimensions

Physics is concerned about the very small, and the very large. In order to represent our observations, we need to know about units. The SI unit for length is the metre (m), time - seconds (s) and mass - kilograms (kg). The L (length), T (time) and M (mass) are fundamental units which all other quantities in physics can be derived.

[ Speed ] = [ L ] / [ T]

The dimension of Speed, is the dimension of Length, divided by the dimension of Time.

[ Volume ] = [ L ]³

The dimension of Volume, is the dimension of Length to the power of 3.

[ Density ] = [ M ] / [ Volume ] = [ M ] / [ L ]³

[ Acceleration ] = [ L ] / [ T ]²

All other quantities can be derived from these 3 fundamentals.

Uncertainty

Any measurement, without any knowledge of its uncertainty, is completely meaningless. The uncertainty is the known error margin. Take for example, two measurements taken are:

25.0cm
24.5cm

Whether it is true that the subject in the first measurement is longer than the second one, or is there a possibility that they might be equal, depends on the error margin or the uncertainty of a measurement. If the uncertainty is +-0.1, we can almost guarantee that the first subject is longer than the second. However, if the uncertainty is +-1cm, we may have to take more measurements.

We can find the uncertainty of our measurements through control tests.

Galileo Proportions Problem

When it comes to measurements in systems, we also need to look at how certain measurements scale with each other. Let's look at an example posed by Galileo Galilei.

Suppose that we have an animal of size S. It has legs with a femur bone of length L. The femur bone has a thickness of D.

It is safe to say that an animal has a femur bone that is proportionate (∝) to its size:
S ∝ L

In this case, it is also safe to say that the mass of an animal is proportional to its size to the power of 3:
M ∝ S³

Therefore:
M ∝ L³

Now notice that in order to support its weight, according to Galileo, its femur bone's cross section must be proportionate to the mass:
M ∝ D²

We can then derive this proportionality:
D ∝ S^3/2

Then he raised a comparison between an elephant and a mouse. An elephant is about 100 times the size of the mouse, so its femur bone would be 100 times longer, which is true. This proportionality suggests that an elephant's femur cross section would also be 1000 times thicker than the mouse, which turned out to be false. The elephant's femur, is only 100 times thicker, which scientists concluded to be nature's protection against buckling.

Derivation of Dimensions

Let us look at how we can derive the dimensions of units. Recall the following formula:

F = ma

[ Force ] = [M] * [ Acceleration ]

[ Acceleration ] = [ Velocity ] / [ T ]

[Velocity ] = [ L ] / [ T ]

Therefore:

[ Acceleration ] = [ L ] / [ T ]²

[ Force ] = [ M ] [ L ] / [ T ]²

Now let's look at Momentum. Momentum is described as:

p = mv

[ Momentum ] = [ M ] * [ Velocity ]

Therefore:

[ Momentum ] = [ M ] * [ L ] / [ T ]

How about Pressure? Pressure is described as Force per unit Area:

[ Pressure ] = [ Force ] / [ Area ]

[ Force ] = [ M ] [ L ] / [ T ]²

[ Area ] = [ L ]²

Therefore:

[ Pressure ] = [ M ] / [ L ] [ T ]²

Finally, Kinetic Energy is defined as:

1/2*mv²

This gives us:

[ Kinetic Energy ] = 1/2 [ M ] [ L ]² / [ T ]²

From here, we can derive Work, which is the application of Force for a certain Distance:

[ Work ] = [ M ] [ L ]² / [ T ]²

We can also derive Power, which is the supply of Energy (which is equal to Work) over Time:

[ Power ] = [ M ] [ L ]² / [ T ]³

Period of Pendulum Problem

Let us now create a function that can tell us the period of a pendulum (time it takes for a full oscillation).

Suppose that we have the following information:
l = Length of Pendulum (L)
m = Mass of Bob (M)
g = Gravitational Acceleration (L/T²)
θ = Angular Amplitude (L)
C = Constant

Our final solution would be in the form:
[T] = C [L]^p [M]^q [L]^r [LT^-2]^s

From here, we can conclude that q is zero because it doesn't appear on the left at all. We can also conclude that s must be -0.5 because T on the left is to the power of 1 only. Since L doesn't appear on the left, then L should be eliminated as well.

The total power of L in this equation is currently p+r+s. In order for L to have the power of 0, we must have:
p+r+s = 0
p+r = -s
p+r = 0.5

Now, in order to find out exactly what p and r is, we'll need to go through experimentation to find out if Angular Amplitude plays a part in the equation. If we mess around with a pendulum, you'll soon discover that the Angular Amplitude does not affect the oscillation at all. Therefore, p = 0.5, r = 0.

Our final expression would be:
[T] = C [L]^0.5 [LT^-2]^-0.5
T = C(l/g)^0.5

If we scale the length by x, the function would react as:
T(scaled) = C(xl/g)^0.5
T(scaled) = (x)^0.5C(xl/g)^0.5
T(scaled) = (x)^0.5T

This is the furthest we can go without experimentation.

Computer Science 17

Noise, in photography, is the inability of your camera's sensor to accurately sample and reproduce a pixel from a given exposure. However, the fact that your picture is still discernible is due to the high Signal-to-Noise Ratio.

Signal-to-Noise Ratio is a measurement of how much of a given piece of information is correct, and how much of it is simply noise. For a typical picture, the SNR is actually extremely high, which means that most of the information is correct.

Noise, like most physical phenomena, is Normally Distributed. In case of a picture, for each pixel, the accuracy of the color representation is normally distributed between 3 axes, the Red, Green and Blue, with the Mean of the distribution being the most accurate color representation of the pixel that a camera can produce. If the color inaccuracy exceeds a certain threshold, we will call it "Noise".

If we measure the Noise of a single image, we measure the number of correctly produced pixels, versus the pixels that are incorrectly produced. We must, of course, define what Noise means in an image. In our example, we'll say that pixels of colors that are below 80% accuracy would be considered Noise. If for every 100 pixels there exists 1 pixel of Noise, we would have the SNR of 100:1. This also means that 1% of total pixels are actually Noise.

In the context of an image, I'll refer to noise as Image Noise. Image Noise of an image is relatively stable. The percentage of pixels that are incorrectly produced for a given exposure is roughly constant and does not fluctuate greatly. When considering Image Noise, an exposure can be 100:1 SNR, the next exposure can be 105:1, and the following can be 95:1.

However, if we look at the noise of each pixel, we measure how far is the color away from the actual. For a 100:1 SNR image, 99% of the pixels are actually sufficiently correctly reproduced, while 1% strayed too far from the Mean.

In the context of individual pixels, I'll refer to the noise as Pixel Noise. Pixel Noise can fluctuate greatly. A pixel can have RGB SNR ranging from 100:1 , or 5:1, or even 1:1 accuracy. It is this fluctuation that causes Pixels to contribute to Image Noise.

In order to fix this, we use Image Averaging. If we average every single pixel across 100 exposures of the same angle, we would have a pixel that is 99% accurate to its actual color. By averaging, no pixels would be 100% accurate this way, but pixels that are below 80% accuracy would become extremely rare. The overall Noise is unchanged, but by evening and spreading out the Pixel Noise across all pixels, we effectively reduced the Image Noise of the image.

Let's go through an example of Image Averaging. You can do this with any series of images taken from the same angle. For our purpose, we'll use frames from a video. Here's the set of frames we have, downsized from 720p:



This is the full image of a single frame. Click to open it in a separate page to further inspect it:



After averaging, this is what we got:



Not bad for a zoomed video at 60x huh?

Let's look at another example, at 30x zoom with some contrasting colors:



Notice that you can also have images where the lighting change slightly. It would all even out perfectly.

(There's been some error with the image hosting. I'll host the images once the server is online)

Here's a single shot in full resolution:



And here's the shot after averaging the pixels across 10 shots:



As you can see, the focal blur is better represented and so on.

This concludes this article. I may write another article on how we can do the averaging, especially across hundreds of images.

Computer Science 16

Greedy Algorithms makes up the answer through Locally Optimal solutions, and therefore is quicker at arriving at an answer. However, it does not guarantee that the final answer is Globally Optimal. The problem, of course, with finding the Globally Optimal is that it is prohibitively expensive to compute. Finding the Global Optimum is Inherently Exponential.

Let us now look at another class of Optimization Problems. One of the most important aspect of Computer Science is the idea of Machine Learning. Superficially, we can define Machine Learning as building programs that learns. However, it suffices to say that almost all programs learn in some form or another.

Machine Learning is a scientific discipline that is concerned with the design and development of algorithms that allow computers to evolve behaviors based on empirical data. We can have a better idea of Machine Learning by looking at the Machine Learning Wikipedia page.

A major focus of Machine Learning is to allow a computer to learn and  recognize complex patterns in order to make decisions based on empirical data. This whole process is known as Inductive Inference. The basic idea is for a computer to observe and record examples from a subset of a statistical phenomena, and generate a model that summarizes some statistical properties of the data in order to predict the future.

There are two distinctive approaches to Machine Learning:
1) Supervised Learning
2) Unsupervised Learning

In Supervised Learning, we associate a Label with each example in a Training Set. If the Label is Discrete, we call it a Classification Problem. For example, we can classify whether a certain transaction on a credit card is done by the owner or not. If the Label is Real, we call it a Regression Problem. When we were doing the Curve Fitting in the previous articles, we were doing Machine Learning, and handling a Regression Problem.

Based on the examples in the Training Set, the goal is create a program in such a way that it is able to predict the answer for other cases before they are explicitly observed. In order to do this, we need to Generalize the statistical properties of the Training Set to be able to make predictions about things we haven't seen.

Let us look at the following Training Set. Suppose that we've obtained some data plotted into axes X and Y, and the Colors (Red and Green) and Shapes (Grass and Leaf) are the Labels of the data.



Upon obtaining a set of data, we must ask ourselves a few questions:
1) Are the Labels accurate?
Labels are subject to machine and human errors.

2) Is the past representative of the future?
If the past results are generally not representative of the future, then it would be extremely difficult to generate a sufficiently accurate model for use.

3) Do we have enough data to generalize?
If the Training Set is relatively small, we should take care not have too much confidence on the predictions.

4) What features are to be extracted?
In order to find out something about the data, we must have the relevant features of the data. For example, in order to find out the composition of haze particles, it is useful to know the types of trees that are being burnt down.

5) How tight should the fit be?
How do we classify the data? Should we fit the data onto the curve, or fit the curve to the data?

It is easy to tell from the Training Data that a certain line y=c separates the labelled Green data from the labelled Red data within the given boundaries.



However, the more difficult part is thinking how should we go about classifying the Shape labels? Here's two of the many possible ways we can classify it:

Do we classify it by grouping the leaves together like this?



This is great because it minimizes the Training Error in the data. However, is this what we really want? How do we know if future data would really fit like this?

Or do we classify it like that, and assume that the leaf with the grass is a labelling error or an outlier? This is a good generalization but there is one Training Error.



That's where we need more data plotted in the vicinity of the stray leaf to find out if it's a Training Error or an Outlier.

Now let's talk about Unsupervised Learning. The big difference here is that we have Training Data, but we don't have Labels. The program gets a bunch of points, but doesn't know what shape or color is it. What can the program learn?

Typically, what the program learns in Unsupervised Learning is the Regularities of Data. From the initial data set, Unsupervised Learning can discover the structures in the data.



The dominant form of Unsupervised Learning is Clustering. What we just did above is to find the Clusters in the data. Clustering is to organize the data into groups whose members are similar to each other. However, how do we describe "similar"? Suppose that the data plots the heights and weights of people, and the x Axis represents weight, while the y Axis represents height. We now have distinct clusters of people.

Walmart in fact used Clustering to find out how to optimally structure their shelves. They found out through Clustering that there is a strong correlation between people who bought beer and people who bought diapers. Therefore they decided to place their beer and diapers departments side by side, which worked surprisingly well.

Similarly, Amazon used Clustering to find people who like similar books based on their profile and buying habits. Netflix uses the same to recommend movies. Biologists use clustering a lot to classify plants and animals. They also use it a lot in genetics to find genes that look and function alike.

What properties do a good Clustering have? It should have:
1) Low Intra-Cluster Dissimilarities
All the points in the same cluster should have as little dissimilarities as possible.

2) High Inter-Cluster Dissimilarities
All points between different clusters should be highly dissimilar.

How do we measure the Intra- and Inter-Cluster Dissimilarities? Well, of course, if you thought of Variance, you are absolutely right!

Intra-Cluster Variance looks at the Variance of the Intra-Cluster Points with respect to the Cluster's Mean.

variance(C) = Sum of, from i=0 to i=len(C), (mean(C)-c_i)²

Inter-Cluster Variance looks at the Variance between the Means of the different Clusters.

Clustering can be formulated as an Optimization Problem. The Optimization Problem of Clustering is: Find the set of clusters, C, such that Intra-Cluster Variance is minimized, and the Inter-Cluster Variance is maximized.

However, such a definition is not enough, because if that's the case, then we can just put each point in its own cluster. The variance would be 0 and Intra-Cluster Variance would be maximized.

What we stated is just the Objective Function. In order to have a proper Optimization Problem, we need to spell out the Constraints. The Constraints can be things like:
1) We can only have a maximum of x clusters
2) There should be a minimum distance between two clusters.

Finding the absolute best solution to this problem is Computationally Prohibitive (it is minimum in the order of n²). Therefore people typically resort to some form of Greedy Algorithm. The two most common Greedy Algorithm used to solve Clustering is:
1) Hierarchical Clustering
2) k-means

Let's go through an example in Hierarchical Clustering. Suppose that we have N items, and we have an accompanying N*M matrix, where M=N, which we can look-up to find out the distance between any two points.

The algorithm is as follows:
1) Assign each item to its own cluster. Therefore, if we have N items, we have N clusters.
2) Find the most similar pair of clusters and merge them.
3) Continue the process until all items are in x number of clusters.

This type of Clustering is also known as Agglomerative Clustering. Step 2 of the algorithm is the most difficult. However, what do we compare if there are two items in each cluster now? What we need to look at is the Linkage Criteria.

There are various Linkage Criteria which we can look at. These are:
1) Single-Linkage
The shortest distance between any two pairs between clusters.

2) Complete-Linkage
The shortest distance between two furthest pairs between clusters. This takes care of the worst case.

3) Average-Linkage
The distance between the means of each cluster.

Of course, like any Greedy Algorithm, there is no Linkage Criteria that's the best. This algorithm is minimum an Order of n² and there is no guarantee that the Cluster is Globally Optimal.

The most important issue in Machine Learning is Feature Selection in order to know what to compare. To Cluster countries together, we need to have something known as a Feature Vector containing all the features that are to be compared. An example of a Feature Vector of a country is a list containing the Population, Coordinates, Size, and so on.

Let's do an example of Hierarchical Clustering. Suppose that we have these coordinates of MRTs:

Dhoby Ghaut,1.298593,103.845909
Buona Vista,1.307412,103.789696
Paya Lebar,1.318089,103.892982
Mountbatten,1.306306,103.882531
Jurong East,1.334308,103.741958
Tampines,1.353092,103.945229
Tanah Merah,1.327257,103.946579
Bishan,1.350772,103.848183
Yio Chu Kang,1.381905,103.844818
Choa Chu Kang,1.385482,103.74424

If we convert the Longitudinal and Latitudinal data into actual distances, we would end up with the following table:

Location	Dhoby Ghaut	Paya Lebar	Tampines	Tanah Merah	Buona Vista	Yio Chu Kang	Bishan	Jurong East	Mountbatten	Choa Chu Kang
Dhoby Ghaut	0	5662	12590	11632	6324	9260	5804	12215	4159	14863
Paya Lebar	5662	0	6989	6043	11540	8885	6163	16880	1750	18148
Tampines	12590	6989	0	2875	18015	11609	10788	22686	8694	22625
Tanah Merah	11632	6043	2875	0	17574	12837	11243	22754	7489	23399
Buona Vista	6324	11540	18015	17574	0	10299	8091	6089	10318	10040
Yio Chu Kang	9260	8885	11609	12837	10299	0	3480	12596	9389	11185
Bishan	5804	6163	10788	11243	8091	3480	0	11946	6244	12179
Jurong East	12215	16880	22686	22754	6089	12596	11946	0	15929	5693
Mountbatten	4159	1750	8694	7489	10318	9389	6244	15929	0	17709
Choa Chu Kang	14863	18148	22625	23399	10040	11185	12179	5693	17709	0

Let's go through this semantically, using the Single-Linkage example:

10 Clusters
Dhoby Ghaut, Paya Lebar, Tampines, Tanah Merah, Buona Vista, Yio Chu Kang, Bishan, Jurong East, Mountbatten, Choa Chu Kang

9 Clusters
Dhoby Ghaut, [Paya Lebar, Mountbatten], Tampines, Tanah Merah, Buona Vista, Yio Chu Kang, Bishan, Jurong East, Choa Chu Kang

8 Clusters
Dhoby Ghaut, [Paya Lebar, Mountbatten], [Tampines, Tanah Merah], Buona Vista, Yio Chu Kang, Bishan, Jurong East, Choa Chu Kang

7 Clusters
Dhoby Ghaut, [Paya Lebar, Mountbatten], [Tampines, Tanah Merah], Buona Vista, [Yio Chu Kang, Bishan], Jurong East, Choa Chu Kang

6 Clusters
[Dhoby Ghaut, Paya Lebar, Mountbatten], [Tampines, Tanah Merah], Buona Vista, [Yio Chu Kang, Bishan], Jurong East, Choa Chu Kang

5 Clusters
[Dhoby Ghaut, Paya Lebar, Mountbatten], [Tampines, Tanah Merah], Buona Vista, [Yio Chu Kang, Bishan], [Jurong East, Choa Chu Kang]

4 Clusters
[Dhoby Ghaut, Paya Lebar, Mountbatten], [Tampines, Tanah Merah], [Yio Chu Kang, Bishan], [Jurong East, Choa Chu Kang, Buona Vista]

3 Clusters
[Dhoby Ghaut, Paya Lebar, Mountbatten, Yio Chu Kang, Bishan], [Tampines, Tanah Merah], [Jurong East, Choa Chu Kang, Buona Vista]

2 Clusters
[Dhoby Ghaut, Paya Lebar, Mountbatten, Yio Chu Kang, Bishan, Jurong East, Choa Chu Kang, Buona Vista], [Tampines, Tanah Merah]

1 Cluster
[Dhoby Ghaut, Paya Lebar, Mountbatten, Yio Chu Kang, Bishan, Jurong East, Choa Chu Kang, Buona Vista, Tampines, Tanah Merah]

As you can tell, at the point where there's 4 Clusters left, it's already looking pretty good. In the next article we will go through the actual implementation of this algorithm!

Computer Science 15

There are many ways we can create models to solve problems. The generalized process is as shown:
1) Start with an experiment in which we gather data.
2) Use computation to find and evaluate a model.
3) Use theory, analysis and computation to derive a consequence of the model.

Optimization Problems involve writing programs to optimize real-life scenarios. Each Optimization Problem consists of an Objective Function and a set of Constraints that has to be satisfied. The Objective Function is used to compare and find the optimal results, while Constraints are rules that must be obeyed (for example, in the case of navigation, Constraints can be in the form of limiting to only pedestrian routes).

Once we come up with these things, we can use computation to attack the problem. A way to solve seemingly new problems is to attempt to map these problems to classic problems that we already know how to solve. This is known as Problem Reduction.

Optimization Problems typically take a long time to solve. Often, there is no computationally efficient way to solve them. Therefore, it is common that we see "Best Effort" results.

A classic Optimization Problem is the Knapsack Problem. The problem is typically discussed in the context of a burglar. One of the problems a burglar must deal with is deciding what to steal. There is usually far more to steal than you can carry away. The Objective Function must optimize the value of what we steal, with the Constraint of the maximum weight or volume we can carry in the knapsack.

Suppose that the burglar has done his research on the house and finds out that these are the objects of interest:

Item	Value	Weight
Desktop	$1500	10Kg
Laptop	$1200	1.5Kg
LED TV	$1000	10Kg
Smartphone	$750	0.25Kg
Ring	$1000	0.01Kg
Piano	$6000	180Kg
HiFi Set	$500	5Kg
PSP	$300	0.5Kg

(Here I'm saying Weight in a general sense, don't tell me it should be Mass instead. I know :p)

One of the easiest solution is to implement a Greedy Algorithm. A Greedy Algorithm is iterative, and at each step, the locally optimal solution is picked. What we need to know is to find out what is considered "locally optimal". It can be the best price/weight ratio, or the highest raw value that can still fit into the bag, or the item with the least weight. There is no guarantee that each variable is the best. Therefore, the biggest problem with the Greedy Algorithm is that we have no single best "locally optimal" that would work most of the time for the 0/1 Knapsack Problem.

The term 0/1 refers to the fact that we can't take half an item. It's either the full item, or we don't take it. Greedy Algorithm, however, is the best for Continuous Knapsack Problem. Take for example, if we break into a smith shop and find barrels of gold, silver and other precious metals. We can then work by filling up with gold first, then silver, then the other metals, and when we finally run out of space, we can take a partial barrel. However, most problems in life are 0/1 Knapsack problems.

Let's create a quick object that allows us to define what an item is:

class Item(Object):
    def __init__(self, name, value, weight):
        self.name = name;
        self.value = float(value);
        self.weight = float(weight);

    def getName(self):
        return self.name

    def getValue(self):
        return self.value

    def getWeight(self):
        return self.weight

    def __str__(self):
        return str(self.name)+", $"+str(self.value)+", "+str(self.weight)+"Kg"

We then instantiate the items based on the above table:

def instantiateItems():
    names=["Desktop","Laptop","LED TV","Smartphone","Ring","Piano","HiFi Set","PSP"]
    values=[1500,1200,1000,750,1000,6000,500,300]
    weights=[10,1.5,10,0.25,0.01,180,5,0.5]
    items=[]
    for i in range(len(names)):
        items.append(Item(names[i],values[i],weights[i]))
    return items

At this point we should have a list of items ready for us to work with. We then come up with the greedy algorithm that does all these. The Greedy Algorithm should accept a list of items, the maximum weight that a person can carry, and the function used for comparison.

Here is the Greedy Algorithm:

def greedy(items, maxWeight, keyFunction):
    itemsSorted = sorted(items, key=keyFunction, reverse=True)
    knapsack=[]
    totalValue=0.0
    totalWeight=0.0
    for item in itemsSorted:
        if totalWeight+item.getWeight()<maxWeight:
           knapsack.append(item)
           totalValue+=item.getValue()
           totalWeight+=item.getWeight()
    return knapsack, totalValue, totalWeight

def value(item):
    return item.getValue()

def weight(item):
    return item.getWeight()

def vwRatio(item):
    return item.getValue()/item.getWeight()

A little explanation here about the formal parameter keyFunction. The Key Function is a function of 1 formal parameter that the sorted() function would use to compare the items in the list. value(), weight() and vwRatio() are examples of Key Functions. Using value() as the keyFunction, for example, would compare whatever is returned by the item.getValue() of all the items.

Let's test the Greedy Algorithm now:

def printResults(knapsack,totalValue,totalWeight):
    for item in knapsack:
        print(item)
    print("Total Value: $"+str(round(totalValue,2)))
    print("Total Weight: "+str(round(totalWeight,2))+"Kg")

def testGreedy():
    items=instantiateItems()
    knapsack,totalValue,totalWeight=greedy(items,22.5,value)
    print("Using value as the key function:")
    printResults(knapsack,totalValue,totalWeight)
    knapsack,totalValue,totalWeight=greedy(items,22.5,weight)
    print()
    print("Using weight as the key function:")
    printResults(knapsack,totalValue,totalWeight)
    knapsack,totalValue,totalWeight=greedy(items,22.5,vwRatio)
    print()
    print("Using vwRatio as the key function:")
    printResults(knapsack,totalValue,totalWeight)

Assuming our burglar has a knapsack that can carry only up to 22.5Kg of things, here's the results:

Using value as the key function:
Desktop, $1500.0, 10.0Kg
Laptop, $1200.0, 1.5Kg
LED TV, $1000.0, 10.0Kg
Ring, $1000.0, 0.01Kg
Smartphone, $750.0, 0.25Kg
PSP, $300.0, 0.5Kg
Total Value: $5750.0
Total Weight: 22.26Kg

Using weight as the key function:
Ring, $1000.0, 0.01Kg
Smartphone, $750.0, 0.25Kg
PSP, $300.0, 0.5Kg
Laptop, $1200.0, 1.5Kg
HiFi Set, $500.0, 5.0Kg
Desktop, $1500.0, 10.0Kg
Total Value: $5250.0
Total Weight: 17.26Kg

Using vwRatio as the key function:
Ring, $1000.0, 0.01Kg
Smartphone, $750.0, 0.25Kg
Laptop, $1200.0, 1.5Kg
PSP, $300.0, 0.5Kg
Desktop, $1500.0, 10.0Kg
LED TV, $1000.0, 10.0Kg
Total Value: $5750.0
Total Weight: 22.26Kg

As you can see, there is no single best Key Function that we can use. In the end it's all up to the situation and we must look at every single one to know what is the best. Since it's called Greedy Algorithm, there is definitely some sort of Order of Complexity associated with it.

The implemented algorithm is divided into two parts: Sorting, then Comparing.

We can assume that the sorting is some variant of the Merge Sort, of which we have:
O(len(items)*log(len(items))) or otherwise O(n log n)

The comparison part is just a for loop that runs a maximum of len(items) times. Therefore we have:
O(len(items)) or otherwise O(n)

Taking the worst case, the Order of Complexity of the algorithm is O(n log n).

However, it is important to note that what we have may not be the absolute optimal we can have.

If we formulate a problem well enough, it would definitely be apparent that it can be solved in a bruteforce way. The biggest trouble is whether it is feasible or not. Fortunately, bruteforce is feasible for small 0/1 Knapsack Problems.

Suppose that we have these 8 items. The maximum combination of items that we can take is:
2⁸ = 256

We can then run a bruteforce algorithm to go through every single possible combination to find out what is the maximum [Sum of Value] where [Sum of Weight] is below maximum weight. In other words:
Objective Function - Maximum [Sum of Value]
Constraint - [Sum of Weight] Below or Equal to Maximum Weight

The Order of Complexity for such an algorithm is:
O(len(items)**2) or O(n²)

Of course, if we're doing just 8 items, we have 256 different combinations. However, what if we have a list of 50 items. That would be 1125899906842624 item combinations. Even if the computer could generate a combination in a nanosecond, it would take 36 years to complete.

Just for kicks, here's the bruteforce algorithm to solve this:

def brange(integer):
    binary=bin(integer)[2:]
    brange=[]
    for i in range(len(binary)):
        if binary[len(binary)-i-1]=='1':
            brange.append(i)
    return brange
   
def bruteForceKnapsack(items,maxWeight):
    combinations=[]
    values=[]
    weights=[]
    maxValue=0.0
    for combination in range(2**len(items)):
        currentWeight=0.0
        currentValue=0.0
        for i in brange(combination):
            currentWeight+=items[i].getWeight()
            currentValue+=items[i].getValue()
        if currentWeight<=maxWeight:
            if maxValue<currentValue:
                combinations=[]
                values=[]
                weights=[]
                maxValue=currentValue               
                combinations.append(combination)
                values.append(currentValue)
                weights.append(currentWeight)
            elif maxValue==currentValue:
                combinations.append(combination)
                values.append(currentValue)
                weights.append(currentWeight)
    return combinations,values,weights

def testBruteForceKnapsack():
    items=instantiateItems()
    combinations,values,weights=bruteForceKnapsack(items,22.5)
    print("For a maxweight of 22.5:")
    for i in range(len(combinations)):
        for j in brange(combinations[i]):
            print(items[j])
        print("Total Value: $"+str(round(values[i],2)))
        print("Total Weight: "+str(round(weights[i],2))+"Kg")
        print()

The brange() function returns a list of 1's in an integer converted to binary. For example, 10 when converted to binary is 1010. It would return 1 and 3 which are the location of the 1's. This algorithm goes through all possible combinations. If there's a better combination than one it currently knows, it would replace all previous combinations. If there's a combination that gives the same value as the one it currently has, it would append it to the list. Here's some output:

For a maxweight of 22.5:
Desktop, $1500.0, 10.0Kg
Laptop, $1200.0, 1.5Kg
LED TV, $1000.0, 10.0Kg
Smartphone, $750.0, 0.25Kg
Ring, $1000.0, 0.01Kg
PSP, $300.0, 0.5Kg
Total Value: $5750.0
Total Weight: 22.26Kg

Here we see that the Greedy Algorithm got lucky. Its Best Effort results were correct! If we work with a different set of items the Greedy Algorithm results may not be as optimal, but the Brute Force it will definitely give you the best, or a set of best results.

Brute Force, of course, is unscalable. In the next article we will look at other ways we can find such results.

Here's graphically how much we can steal from the house:



As you can see, optimally 28Kg is minimum to steal mostly everything. You'll have to add more than 150Kg to your maximum load in order to start steal more.

Computer Science 14

Before believing the result of any simulation, we have to have confidence that our conceptual model is correct, and that we have correctly implemented the model.

Before, we looked at Statistical Tests. These tests give us Statistical conclusions which tell us things about the simulation, about how it converges to a certain result. However, we still have to perform Sanity Checks.

One of the ways we can do this is to test the results against reality. For example, in the previous article, the value of Pi obtained from the simulation can be tested with the formula to know if we got the right answer or not.

Let us look at an interplay between Physical Reality, Theoretical Models and Computational Models. The Physical Reality is based on the Physical Situation. An example of a Physical Situation can be the price of Big Mac across countries (which we know as the Big Mac Index), or a stock market, or viruses in a body. We use a Theoretical Model to have some deeper insights into it, and when it gets too complicated or is unable to show us everything, we use a Computational Model.

At times when we conduct an actual experiment, we may not get the results that match exactly with the theory. To work with experimental errors, we have to assume that there is some sort of random perturbation applied to the results. Gaussian's work suggests that we can model experimental errors as Normally Distributed.

The Hooke's Law, F=kx, was described by Robert Hooke in 1678. In plain English, the Hooke's Law states that the Force (F) stored in a spring is linearly related to the Distance (x) the spring has either been compressed or stretched multiplied by the Stiffness (k) of the spring. A stiff spring has a large k value. This law holds for a wide variety of materials and systems. However, of course, it does not hold for a force that causes the spring to exceed its elastic limit.

We can find the k value of a spring by suspending the spring and hanging known masses on it. Measuring the distance for each mass. Suppose that a 1K
Using the rules:
F = kx
F = ma
ma = kx
k = ma/x

We derive:
k
= 1Kg * 9.81m/s² / 0.2m
= 49.05N/m

We can do a sanity check to make sure this is correct (using 0.4m as the distance, we should evaluate to 2Kg):
49.05 = m * 9.81m/s²  / 0.4m
49.05 / (9.81m/s²/0.4) = 2Kg

What we are assuming here is that 49.05N/m is the correct k value with just one experiment. However, remember that experimental errors are normally distributed. If we actually put a 2Kg weight there, we may or may not get a stretching of exactly 0.4m. People would typically repeat the same experiment many times with different weights in order to get enough data to compute an accurate value of k.

Suppose that we have the following data recorded by a fellow scientist:

0.0 0.0
0.05 0.00968
0.1 0.02087
0.15 0.02871
0.2 0.03849
0.25 0.05371
0.3 0.05917
0.35 0.07066
0.4 0.08041
0.45 0.08863
0.5 0.10219
0.55 0.10251
0.6 0.09903
0.65 0.10078
0.7 0.10321
0.75 0.09605
0.8 0.09730
0.85 0.09986
0.9 0.09876
0.95 0.09825
1.0 0.09986

The data comes in a file called spring_data.txt describing weights (in kg) and stretch (in m). In order to manipulate this, we can either manually copy the data into our program, or import and parse it from within Python. In this example, we are going to import and parse the data:

The following code imports the data, and then returns two lists:

def importData(fileName):
    masses=[]
    distances=[]
    file=open(fileName,'r')
    for line in file:
        m,d=line.split()
        masses.append(float(m))
        distances.append(float(d))
    file.close()
    return masses,distances

masses,distances=importData("spring_data.txt")
print(masses)
print(distances)

Now that we have the two lists, we can then plot the data using:

import pylab

pylab.plot(masses,distances,"o")
pylab.xlabel("Mass (kg)")
pylab.ylabel("Distance (m)")
pylab.title("Spring Experiment")
pylab.show()

As you can see, the distances are somewhat linear at first (with some level of normally distributed error), then it evens off at the end:



In Numpy, there is a type called the "array". This array has nothing to do with the array we have in Java. The array here can be considered as more of a matrix. Changing data to an array allows:
1) Mathematical operations to be performed on it (i.e. array*3 will multiply every element by 3)
2) Cross-products between two arrays (array1*array2)

Lets look at some application of arrays:

import pylab

forces = pylab.array(masses)*9.81

pylab.plot(forces,distances,"o")
pylab.xlabel("Force (N)")
pylab.ylabel("Distance (m)")
pylab.title("Spring Experiment")
pylab.show()

In this example, the masses list is changed into an array and then multiplied by 9.81. What we get is a list of Force (N):



From here, we can calculate k. But before that,we must know whether our data is sensible. In our Theoretical Model, the data should fall on a line F=kx. If we can draw a good line here, we can find k very easily.

If we only have two points, we can easily join them to create a line. If we have a bunch of points scattered around, we must find the line that is the closest to all the points. This is known as Fitting. In order to find the line, we need to have a measure to how good is the fit. We need an Objective Function to compare fits. The standard Objective Function is the Least Squares Fit.

Suppose we have two lists, a list of predictions according to the line, and a list of observations. The Least Squares Fit is described as below:

def leastSquaresFit(observed,predicted):
    sum=0.0
    for i in range(len(observed)):
        sum+=(predicted[i]-observed[i])**2
    return sum

For comparing two Line of Best Fits (or Polynomial of Best Fits), the lower the number, the better.

A good model is created when, for each independent variable Mass (kg), we can predict a good dependent variable Displacement (m).

For creation of the Line of Best fit of any graph, we use Linear Regression. Pylab does Linear Regression using the function polyfit(observedX,observedY,degreeOfPolynomial). Recall that a line is y=mx+c. Polyfit would return the values of m and c for 1 degree fitting:

import pylab

forcesArray = pylab.array(masses)*9.81
distancesArray = pylab.array(distances)

m,c = pylab.polyfit(forcesArray,distancesArray,1)
distancesFitArray = m*forcesArray+c

pylab.plot(forcesArray,distancesArray,"o")
pylab.plot(forcesArray,distancesFitArray)
pylab.xlabel("Force (N)")
pylab.ylabel("Distance (m)")
pylab.title("Spring Experiment, k="+str(1/m))
pylab.show()

We have the following output:



We can compute for k through this reasoning:
y=mx+c
Distance = m * Force + c
Force =  1/m * Distance + c
k = 1/m

So is the k value correctly shown? If we look at the picture, we'll realize that the line is not a very good fit. Is it better to use a cubic fit to model our data? Let's look at the example below:

import pylab

forcesArray = pylab.array(masses)*9.81
distancesArray = pylab.array(distances)

m,c = pylab.polyfit(forcesArray,distancesArray,1)
distancesFitArray = m*forcesArray+c
a,b,c,d = pylab.polyfit(forcesArray,distancesArray,3)
distancesFitCubicArray = a*forcesArray**3 + b*forcesArray**2 + c*forcesArray**1 + d

pylab.plot(forcesArray,distancesArray,"o")
pylab.plot(forcesArray,distancesFitArray)
pylab.plot(forcesArray,distancesFitCubicArray)
pylab.legend(["Observed","Linear Fit","Cubic Fit"])
pylab.xlabel("Force (N)")
pylab.ylabel("Distance (m)")
pylab.title("Spring Experiment, k="+str(1/m))
pylab.show()

It seems as though the Cubic Fit is much better than the Linear Fit in describing the spring?



The useful thing about a model is its ability to predict data. Let's look at how the Cubic Fit would predict our data. We want to see, for example, what would happen if an approximately 20N force is placed on the spring. We do this by adding another point in the Forces array:

import pylab

forcesArray = pylab.array(masses)*9.81
forcesArrayExtended = pylab.array(masses+[2])*9.81
distancesArray = pylab.array(distances)

m,c = pylab.polyfit(forcesArray,distancesArray,1)
distancesFitArray = m*forcesArray+c
a,b,c,d = pylab.polyfit(forcesArray,distancesArray,3)
distancesFitCubicArray = a*forcesArrayExtended**3 + b*forcesArrayExtended**2 + c*forcesArrayExtended**1 + d

pylab.plot(forcesArray,distancesArray,"o")
pylab.plot(forcesArray,distancesFitArray)
pylab.plot(forcesArrayExtended,distancesFitCubicArray)
pylab.legend(["Observed","Linear Fit","Cubic Fit"])
pylab.xlabel("Force (N)")
pylab.ylabel("Distance (m)")
pylab.title("Spring Experiment, k="+str(1/m))
pylab.show()

Now observe the output:



The Cubic Fit seems to suggest that if we place 20N of force on the spring, the spring would retract back higher than it originally was. That happens maybe in another universe with different physics but that doesn't happen on Earth.

Remember that Hooke's Law is a first order linear approximation that will fail as the spring approaches its limit of elasticity. Just because the data is not plotted into a straight line does not mean that we should attempt to fit an arbitrary curve onto it, because any polynomial with a high enough degree can fit into any data. What we want, instead, is to look at the portions where Hooke's Law remains valid - that is, when the spring hasn't hit its limit of elasticity.

Let us repeat the experiment, discarding the last 10 results:

import pylab

forcesArray = pylab.array(masses)*9.81
forcesArrayLimited = pylab.array(masses[:-10])*9.81
#forcesArrayExtended = pylab.array(masses+[2])*9.81
distancesArray = pylab.array(distances)
distancesArrayLimited = pylab.array(distances[:-10])

m,c = pylab.polyfit(forcesArrayLimited,distancesArrayLimited,1)
distancesFitArray = m*forcesArray+c
#a,b,c,d = pylab.polyfit(forcesArray,distancesArray,3)
#distancesFitCubicArray = a*forcesArrayExtended**3 + b*forcesArrayExtended**2 + c*forcesArrayExtended**1 + d

pylab.plot(forcesArray,distancesArray,"o")
pylab.plot(forcesArray,distancesFitArray)
#pylab.plot(forcesArray,distancesFitCubicArray)
pylab.legend(["Observed","Linear Fit","Cubic Fit"])
pylab.xlabel("Force (N)")
pylab.ylabel("Distance (m)")
pylab.title("Spring Experiment, k="+str(1/m))
pylab.show()

Now notice that we've obtained a much more believable fit.



Now notice that we've obtained a better approximation for the k value of the spring, which is originally rated at 49.05 by the manufacturer.

To know how good a fit is to a particular set of results, we use the Coefficient of Determination. The Coefficient of Determination is:
r²=1-(Estimated Error / Variance of Measured Data).

Computationally, it is expressed as such:

def rSquare(predicted,measured):
    estimatedError = ((predicted-measured)**2).sum()
    measuredMean = measured.sum()/float(len(measured))
    measuredVariance = (measuredMean-measured**2).sum()
    return 1-estimatedError/measuredVariance

Thanks to Numpy's array, we performed array operations in one line each. Do you see how useful arrays are!?

We can compare the goodness of each fit by comparing it to the measured data. The closer it is to 1, the better. A CoD of 1 means that the fit maps directly onto the measured data.

Computer Science 13

How do we construct computational models that will help us understand the real world?

A Gaussian Distribution, or commonly known as Normal Distribution or the Bell Curve, can be described with just its Mean and Standard Deviation.

Whenever possible, we would model distributions as Gaussian because they are extremely easy to deal with. We have things like the Empirical Rule to help us draw conclusions. However, if a certain distribution is not Gaussian and we treat it that way, the results could be extremely misleading.

Take for example the results of rolling a die. The outcome of die-rolling is not normally distributed. In the case of a fair-die, each face is equally likely to appear. Plotting the results on a histogram would not yield a Bell Curve.

In a fair lottery system, the probability of anything coming up is the same. Most people would study past results in order to predict what comes up next. However, illustrating the Gambler's Fallacy again, lottery has no memory and is thus not based on past results.

In case of the Lottery and Die Rolling, we have a Uniform Distribution. Uniform Distributions describe distributions where each result is equally plausible. The only thing required to describe a Uniform Distribution is the Range of the results.

The Uniform Distribution does not occur much in nature. It is typically observed in fair chance games deviced by humans. Uniform Distribution is not as useful in modelling complex systems as the Gaussian Distribution.

Another type of distribution that occurs quite frequently is the Exponential Distribution. Exponential Distribution is the only Continuous Distribution that is memory-less. The concentration of drug in the bloodstream can be modelled with Exponential Distribution.

Assume that each time-step, each molecule has a probability p of being cleared by the body. The system is memory-less, in the sense that the probability of each molecule being cleared is not dependent on what happened to the previous molecules in the previous steps.

If the probability of being cleared at each step is p, then at time T=1, the probability of the molecule still being in the body is 1-p. Since it is memory-less, the probability of it being cleared is still 1-p. Therefore, the probability of it still being in the body is (1-p)².

We can see that the probability of the molecule still being in the body at the particular time is (1-p)^t

Suppose there are 1000 molecules, and at each time-step, each molecule has a chance of 0.01 to be cleared by the body. We would end up with the following graphs:



Of course, what we are seeing is the mathematical graphing of what we would expect. This graph is achieved through the following code:

import pylab

def mathematicalMolecules(n, clearProbability, steps):
    moleculesRemaining=[n]
    for i in range(steps):
        moleculesRemaining.append(n*((1-clearProbability)**i))
    return moleculesRemaining

moleculeList = mathematicalMolecules(1000,0.01,1000)
pylab.plot(moleculeList)
pylab.xlabel("Time (T)")
pylab.ylabel("Molecules")
pylab.title("Exponential Decay of Molecules")
pylab.show()

Now, notice that the rate of clearing of molecules decrease in a logarithmic way. The less molecules there are, the less molecules are going to be lost. If we plot the graph with semilogy, we can see that it is indeed logarithmic:



We can actually create a Monty Carlo simulation to mimick the process described, instead of simply using mathematical formulas to do it. We make use of the following code to implement this:

import random

def simulationMolecules(n, clearProbability, steps):
    moleculesRemaining=[n]
    for i in range(steps):
        for j in range(n):
            if random.random()<clearProbability:
                n-=1
        moleculesRemaining.append(n)
    return moleculesRemaining

moleculeList = simulationMolecules(1000,0.01,1000)
pylab.plot(moleculeList)
pylab.xlabel("Time (T)")
pylab.ylabel("Molecules")
pylab.title("Exponential Decay of Molecules (Simulation)")
pylab.show()

Notice that the output is not as smooth as we want it to be, but it matches the mathematical diagram almost exactly:



If we look at this graph in a logarithmic way, this is what we'll get:



The biggest difference is that our simulation did not allow for fractions of a molecule, which is more true to what we are looking for.

Here's a direct mapping of the two graphs:





The Mathematical Model is also known as the Analytic Model. Here, we've illustrated the difference between Analytical and Simulation for Exponential Decay. There is no right or wrong answer over here. When looking at the models, we look at its Fidelity and Credibility. Fidelity refers to whether the outcome is accurate. Credibility is whether the outcome can be believed.

Then there is Utility. Which model can be applied to what kind of questions? There is an additional Utility offered by the Simulation Model, because we can easily change the model to be slightly different in ways that are usually harder for an Analytic Model.

For example, if we wish to calculate the rate of which the body clears bacteria, but also factor in things like the rate at which the bacteria can regenerate themselves, it is easier to model with a Simulation rather than an Analytic Model.

Let us look at the Monty Hall problem. In the Monty Hall problem, we start with 3 doors. Behind one of the doors, there is a great prize. Behind each of the other doors, there is a booby prize. A person is asked to choose one of the doors. The host then reveals one of the other two doors with the booby prize. There is then a question posed to the person: Do you want to switch?

The question is: Which is better? Does it matter whether we switch?

Let us look at this first in an Analytical perspective:
The chance that the prize lies in your door is 1/3
The chance that the prize lies in one of the other 2 doors is 2/3

After choosing the first door, one of the two other doors is eliminated. Now, what you have is:
The chance that the prize lies in your door is 1/3
The chance that the prize lies in the other (2 - 1) door is 2/3

Therefore, making the switch doubles your chance from 1/3 to 2/3.

The following maps the problem out in a computational way:

import random
import pylab

def montyHall(n):
    won=0
    for i in range(n):
        correctDoor = random.randint(1,3)
        chosenDoor = random.randint(1,3)
        if correctDoor==chosenDoor:
            won+=1
    return won

def montyHallSwitch(n):
    won=0
    for i in range(n):
        correctDoor = random.randint(1,3)
        chosenDoor = random.randint(1,3)
        if correctDoor!=chosenDoor:
            won+=1
    return won

def plotMontyHall(trialsExp):
    participants=[]
    montyHallList=[]
    for i in range(trialsExp):
        montyHallList.append(montyHall(2**i))
        participants.append(2**i)
    montyHallSwitchList=[]
    for i in range(trialsExp):
        montyHallSwitchList.append(montyHallSwitch(2**i))
    pylab.plot(participants,montyHallList)
    pylab.plot(participants,montyHallSwitchList)
    pylab.xlabel("Participants")
    pylab.ylabel("Winners")
    pylab.title("Monty Hall: Stand vs Switch")
    pylab.legend(["Stand","Switch"])
    pylab.show()

plotMontyHall(15)

This is based upon the following:
1) If guessed door is prize door, stand will give a win
2) If guessed door is NOT prize door, switch will give a win

This is the resulting graph. As we can see, the winning fraction is exactly twice for Switch. Also, 1/3 of Stand participants and 2/3 of Switch participants are winners.:



We can use randomized algorithms to solve problems in which randomness plays no roles. This is surprising, but it is an extremely important tool.

Let us consider Pi (π). Pi is a constant that people have known since the 18th century. The earliest estimate of Pi was 4*(8/9)²=3.16 by the Egyptians in 1650 B.C. In 650 B.C., when Solomon was creating his Basin, he made use of Pi, which was roughly 3. The best estimate of Pi in ancient time was by Achimedes. He didn't give the value of Pi. Instead, he built a polygon of 96 sides, in which he concluded Pi to be between 223/71 to 22/7.

Later on, Buffon and Laplace came up with a Stochastic Simulation to solve for Pi. This is the Buffon and Laplace solution to finding Pi:



Suppose that we have a circle of radius 1 inside a square of sides 2.

We know for sure that the area of the circle is πr², and the area of the square is r²

Buffon and Laplace proposed to randomly drop needles into such a set-up on the ground, then use the number of needles in the circle to estimate the area of the circle and the total number of needles to estimate the area of the square.

We can then use the following to obtain π:
Area of Circle / Area of Square
= πr²/(2r)²
= πr²/4r²
= π/4

π = 4 (Area of Circle / Area of Square)

Therefore:

π = 4 (Needles in Circle / Needles in Square)

Because Buffon and Laplace didn't have the proper physical setup to conduct such a test, we can help him through the use of a simulation:

import math
import random

def estimatePi(n):
    circle=0
    for i in range(n):
        x = random.random()
        y = random.random()
        if math.hypot(x-0.5,y-0.5)<0.5:
            circle+=1
    return 4*circle/float(n)

print (estimatePi(50000000))

In my example, I used a circle of radius 0.5 instead. The formula is the same. We drop 50 million pins into the setup, counted the pins in the circle, and used the formula to give us the answer. The pins were just a means of estimating area.





As you can see, as the amount of pins used increased, the results were closer to Pi and its variation decreased. Here is the code that generated these graphs:

import pylab

def standardDeviation(X):
    mean=sum(X)/float(len(X))
    total=0.0
    for x in X:
        total+=(x-mean)**2
    return (total/len(X))**0.5

def mean(X):
    return sum(X)/float(len(X))

def plotEstimatePi(pinsExp, trials):
    pis=[]
    pins=[]
    stdDev=[]
    for i in range(pinsExp):
        pins.append(trials*2**i)
        currentPis=[]
        for j in range(trials):
            currentPis.append(estimatePi(2**i))
        pis.append(mean(currentPis))
        stdDev.append(standardDeviation(currentPis))
    pylab.plot(pins,pis,"o")
    pylab.xlabel("Pins Used")
    pylab.ylabel("Estimated Pi")
    pylab.title("Pi Estimation, Final Value = "+str(pis[len(pis)-1]))
    pylab.semilogx()
    pylab.figure()
    pylab.plot(pins,stdDev)
    pylab.xlabel("Pins Used")
    pylab.ylabel("Standard Deviation")
    pylab.title("Std Dev, Final Value = "+str(stdDev[len(stdDev)-1]))
    pylab.semilogx()
    pylab.show()
  
plotEstimatePi(15,100)

Note that we can actually come up with more accurate representations of Pi through other means, like binary search. However, this example illustrates the usefulness of randomness in solving non-random problems. This method is extremely useful when solving for integration, because we know that the integration of a function is the area under the graph.