text:
# Thread: Where does the money come from?

1. NNY
Member
Join Date
Mar 2018
Posts
49

## Where does the money come from?

One thing that hasn't been mentioned yet is, where does all the money come from?

The State creates very little of the money that exists in an economy.

I'm not an economist but my understanding is, the Blue line is M0, which is money created by the federal reserve. And everything above that, M1 and M3, is created by the banks from deposits.

The way deposits are turned into more money is by Money Multiplying :
But, this article from The Australian explains, what is supposed to happen, doesn't happen.

So where is all the money that is created by the banks?
Well, it's not in my bank account, that's for sure.

Whoever controls the banks most probably has control over all the money.

I remember a conversation with a lawyer, who worked for some social/human rights group, dealing mostly with cases for indigenous people in Sydney. He told me, if you want to become rich by way of lawyering, become a lawyer for a bank.

NNY, I'm not an economist either. Banks keep money circulating around the economy. If you have $1,000 under your mattress that money is out of the economy. If you have$1,000 in the bank and the bank lends $900 to someone else, then effectively you have$1K and someone else has $900. They pay for something, the payee puts that money in their bank account, the bank lends it out again. It's the same money but the bank keeps it circulating. A reserve bank or central bank can issue genuinely new money. Obviously we know what happens if a central bank issues too much new money. You can do thought experiments by imagining you decide to be the "Smith Street Bank" for example, and your neighbours give you money to look after for them, which you lend out. It's fairly easy to understand that way. 3. NNY Member Join Date Mar 2018 Posts 49 ## Re: Why income/wealth equality is getting worse Hi 142857, The maths for Money Multiplication should be as follows; Deposit$1000, $100 goes to reserve,$900 lent out
that $900 is deposited,$90 goes to reserve, $810 lent out that$810 is deposited, $81 goes to reserve,$729 lent out
that $729 is deposited,$73 goes to reserve, $656 lent out and so on. Effectively, if somebody deposits$1000, the banks can then lend $10,000 based on that$1,000 deposit. And that $10,000 then becomes part of the money circulating in the economy. But The Australia article says that that isn't what is happening. A far greater amount of money, far greater than$10,000, is being lent out based on the original $1,000. Last edited by NNY; 11th April 2018 at 12:47 AM. 4. ## Re: Why income/wealth equality is getting worse NNY said Hi 142857, The maths for Money Multiplication should be as follows; Deposit$1000, $100 goes to reserve,$900 lent out
that $900 is deposited,$90 goes to reserve, $810 lent out that$810 is deposited, $81 goes to reserve,$729 lent out
that $729 is deposited,$73 goes to reserve, $656 lent out and so on. Effectively, if somebody deposits$1000, the banks can then lend $10,000 based on that$1,000 deposit. And that $10,000 then becomes part of the money circulating in the economy. But The Australia article says that that isn't what is happening. A far greater amount of money, far greater than$10,000, is being lent out based on the original $1,000. Yes, good point. I looked at your other links but not at the Australian article until just now. The thing about banking these days is that "household deposits" (which are subject to reserve requirements in some countries) make up only a small portion of what the banks lend out. Banks have far more creative ways of raising funds for lending now. So the multiplier effect can be a lot more than$10K on most of the money that banks lend out.

Remember when banks started putting fees on savings accounts? This was because most savings accounts cost a bank more to administer than the bank can earn by lending that money out, and banks are actively discouraging accounts with small balances. That is a big shift that has happened within the past couple of generations.

This article is a little bit easier to follow, I thought, than the Australian article: http://positivemoney.org/how-money-w...about-banking/

My professional background is mostly on the systems side of lending and accounting at banks, more in the nuts and bolts side than the conceptual side which is what these articles are dealing with. So I'm likely wrong on some of it.

One thing that people get mixed up is the idea that banks are allowed to "create" money to lend out. This leads to the misconception that banks pluck money out of thin air. In reality every dollar that a bank lends out has to come from somewhere. It's the financial system that generates money by recirculating deposits. I think economists tend to do a lousy job of explaining that side of it to the rest of us.

NNY, how do you think all of this ties into the increase in wealth inequality? I can think of one way off the top of my head which I'll expand on after I get some sleep.

5. ## Re: Why income/wealth equality is getting worse

Money multiplication and its contribution to increasing wealth inequality:

Back when we were kids, people put money in a savings account. Banks lent this money out and added a margin of about 3% to the interest rate. The bank employed tellers, accountants, security guards etc. And the traditional money multiplier effect applied.

Now... most people don't save much, not compared to the vast amount of money that banks lend out for things like mortgages these days. But we have forced savings, and people who have been accumulating super for 2 or 3 decades often have hundreds of thousands of dollars in their super accounts. The money in your super goes back into the economy, used to fund bank lending among other things. And there's no reserve to worry about.

The people who run those super funds, who invest your super, are taking a hefty cut and joining the ranks of the wealthy. They don't call Macquarie Bank, for example, the "Millionaire Factory" for nothing.

6. NNY
Member
Join Date
Mar 2018
Posts
49

## Re: Why income/wealth equality is getting worse

142857 said
NNY, how do you think all of this ties into the increase in wealth inequality? I can think of one way off the top of my head which I'll expand on after I get some sleep.
With The Australian newspaper article and the Wikipedia article on Money Multiplier, there is the assumption that M0, being money created by the reserve bank, is the only money that affects M1 and M3. What the 2 articles don't take into account is massive investments from overseas.

In the link, $3.2 trillion was invested in Australia from overseas in 2016. If that money was initially deposited into a bank, then step one of the Money Multiplier would allow banks to 'create'$320 billion on that deposit that could be used for loans, with a potential of a maximum $32 trillion being created. So who is getting a share of that$320 billion to $32 trillion in the form of loans? I'm sure the reserve bank keeps an eye on all the money from overseas being deposited in Australian bank accounts, in order to keep inflation under control and stop money laundering etc. But, I'm also sure I will never find out who is receiving massive loans that were generated from that$3.2 trillion from overseas investors.

I have met one Sydney lawyer who charges $6,000 per hour as an arbitrator for large companies when they have conflicting opinions about contracts they've entered into. It seems a lot, but if he can solve a problem in 5 hours ($30,000) it is still much much cheaper for the companies to use him than going to court for several years. (He doesn't earn that much money constantly everyday though)

So, with massive amounts of overseas investment and the massive amounts of money 'created' from those investments via the money multiplier, I can see why lawyers who work for banks can earn lots of money.

Finally, I can't say the opinions I have written are true. They are just the only plausible reasons I can think of as to why M1 and M3 in the chart are so high compared to M0.

7. AFA Member
Join Date
Oct 2012
Posts
4,788

NNY said

In the link, $3.2 trillion was invested in Australia from overseas in 2016. The 3.2 trillion is not per year but rather like a capital deposit which fluctuates up or down from year to year, with withdrawals or deposits. I'm not quite sure how you're understanding it. Last edited by wadaye; 12th April 2018 at 12:49 AM. 8. NNY Member Join Date Mar 2018 Posts 49 ## Re: Why income/wealth equality is getting worse wadaye said The 3.2 trillion is not per year but rather like a capital deposit which fluctuates up or down from year to year, with withdrawals or deposits. I'm not quite sure how you're understanding it. In the article from the Australia newspaper, it stated the M1 and M3 are far far greater than they should be compared to the M0. So where does this disparity come from? My guess is overseas investments coming into the country which are initially placed into bank deposit accounts. Looking at the stats of the link. 2014 had$2.8 trillion in overseas investments. Which could be used as 280 billion in loans by way of the first step in the Money Multiplier and potentially 28 trillion using the full effect of the Money Multiplier.

2015 had $3 trillion in overseas investments. Which could be used as 300 billion in loans by way of the first step in the Money Multiplier and potentially 30 trillion using the full effect of the Money Multiplier. 2016 had$3.2 trillion in overseas investments. Which could be used as 320 billion in loans by way of the first step in the Money Multiplier and potentially 32 trillion using the full effect of the Money Multiplier.

I'm sure the Reserve Bank of Australia limits how the banks use the Money Multiplier though. If rules and regulations weren't put into effect, inflation and other abuses would be out of control.
Last edited by NNY; 12th April 2018 at 07:55 AM.

9. ## Re: Why income/wealth equality is getting worse

NNY said
I'm sure the Reserve Bank of Australia limits how the banks use the Money Multiplier though. If rules and regulations weren't put into effect, inflation and other abuses would be out of control.
The money multiplication effect is not really a "tool" used by the banks to make money. It's more a by-product of the financial system. Banks make money by raising money at a lower rate and lending it at a higher rate, although that is the easy part of the equation and if that was all there was to it then you or I could start a bank and make lots of money. Risk management is really what banks traditionally need to be good at, although if financial institutions can make money off lending margins while getting someone else to wear the risk.... that's even better (that was what the sub prime mortgage disaster was about).

That $3.2 trillion foreign investment is not subject to fractional reserves, which is the traditional control that the reserve bank has over the money multiplication effect. Overseas investors are buying land, buying equipment, paying wages and all of that generates economic activity and more money obviously. We are part of a global economy, and Australians invest large sums of money overseas as well. 10. ## Re: Why income/wealth equality is getting worse 142857 said The money multiplication effect is not really a "tool" used by the banks to make money. It's more a by-product of the financial system. Banks make money by raising money at a lower rate and lending it at a higher rate, although that is the easy part of the equation and if that was all there was to it then you or I could start a bank and make lots of money. Risk management is really what banks traditionally need to be good at, although if financial institutions can make money off lending margins while getting someone else to wear the risk.... that's even better (that was what the sub prime mortgage disaster was about). That$3.2 trillion foreign investment is not subject to fractional reserves, which is the traditional control that the reserve bank has over the money multiplication effect. Overseas investors are buying land, buying equipment, paying wages and all of that generates economic activity and more money obviously. We are part of a global economy, and Australians invest large sums of money overseas as well.
Money isn't created by borrowing at one rate and lending at a higher one. Because if I deposit money in the bank, and it's then lent to you, you have to get that money from somewhere. It's a zero sum game, if I have a dollar more, then someone has a dollar less.

The only way that you can get more money is by depriving it to someone else.

At least that's my simplistic understanding.

#### Posting Permissions

• You may not post new threads
• You may not post replies
• You may not post attachments
• You may not edit your posts
•


text:
## The Annals of Probability

### The Asymptotic Distribution of Trimmed Sums

#### Abstract

Let $X_{1,n} \leq \cdots \leq X_{n,n}$ be the order statistics of $n$ independent and identically distributed random variables and let $m_n$ and $k_n$ be positive integers such that $m_n \rightarrow \infty, k_n \rightarrow \infty, m_n/n \rightarrow 0$ and $k_n/n \rightarrow 0$ as $n \rightarrow \infty$. We find a necessary and sufficient condition for the existence of normalizing and centering constants $A_n > 0$ and $B_n$ such that the sequence $T_n = A^{-1}_n \big\{\sum_{i=m_n+1}^{n-k_n} X_{i,n} - B_n\big\}$ is stochastically compact and completely describe the possible subsequential limiting distributions. We also give a necessary and sufficient condition for the existence of $A_n$ and $B_n$ such that $T_n$ is asymptotically normal. A variant of Stigler's theorem when $m_n/n \rightarrow \alpha$ and $k_n/n \rightarrow 1 - \beta$, where $0 < \alpha < \beta < 1$, is also obtained as a by-product.

#### Article information

Source
Ann. Probab., Volume 16, Number 2 (1988), 672-699.

Dates
First available in Project Euclid: 19 April 2007

https://projecteuclid.org/euclid.aop/1176991780

Digital Object Identifier
doi:10.1214/aop/1176991780

Mathematical Reviews number (MathSciNet)
MR929070

Zentralblatt MATH identifier
0647.62030

JSTOR

format: plain


text:
# Blog

## Farewell to astronomy

Recently I took the decision to shift my life onto a new track and accepted a job offer outside of astrophysics. Given that leaving academia has long been a frequent source of conversation among the postdocs and PhD students I interact with on a daily basis I thought that --- inspired by Marcel Haas' article titled Leaving the field, becoming an extronomer --- I'd write a few words about why I'm leaving, and where I'm going.

Jessica Kirkpatrick has written a lot about the difference between academia and industry (Astronomer to data scientist, Astronomy vs. data science) and I'd encourage anybody who is interested to read those articles. They reflect my feelings about academia more coherently than I can myself. Therefore, I'm only going to cover my personal circumstances and say that as much as I love doing science for a living, I do feel that I'll be just as fulfilled elsewhere.

### Where have I been?

I have been incredibly fortunate throughout my career. During my PhD at Durham University and post-doctoral positions at Leiden and Chicago I have been surrounded by wonderful and caring people who, over the course of each contract, have become close friends. Since moving to Chicago I have married a wonderful woman, we are lucky enough to live close to her family (this is by no means to be expected for those of us on short-term academic contracts!), and I have found a group of friends who make me incredibly happy to be around.

The process of moving institutions every few years is both an up- and down- side of early-career academia. It's great to meet new people, and to experience new countries. At the same time, it is exhausting. Each time I move, the process of meeting a whole new circle of friends and stocking another apartment with Ikea furniture feels less like fun and more like a horrible chore, especially with the knowledge that my time at the new institution has an expiration date just a few years in the future.

If I try to project my career forwards and objectively place myself next to others who would be competing with me for jobs, I do feel like I would likely be able to find some sort of faculty position. However, it is unlikely that I'd be able to land something at a prestigious research institution, or in a city that I'd love as much as I do Chicago. I'd probably either end up stuck in a holding-cycle where I take more postdoctoral positions in new cities until finally settling, or I would accept a faculty position in a city and department where the things I have grown to love about both the big-city and big-department would be largely unavailable to my family and me. Either way, continuing an astronomy career would, for me, involve more years of uncertainty and big life changes with no obvious path to resolution.

It is a decision that is personal to everybody, but, for me, when I weigh my academic career against my personal happiness I know it's time for me to make some decisions that prioritize my family and my non-academic life above chasing a tenure track position.

### Where am I going?

A while ago I filled out a LinkedIn profile with all of the computer and data-analysis experience that I gained doing astronomy. As a result, I get the occasional call or email from a recruiter asking if I'd be interested in opportunity X or job Y. Up until now I have always politely brushed these offers off because nothing appealed to me more than spending more time doing science, even though I already knew that in the long-term academia would not be for me.

However, a few weeks ago I received one of these emails from a recruiter on behalf of a company called Narrative Science, with the description that Narrative Science was born out of Northwestern University and "Our artificial intelligence platform automatically analyzes data and produces narratives that are contextually relevant, actionable and tailored to any audience." Basically, they turn tables full of data into stories written in the English language. Narrative Science made a bit of a splash a few years ago with software that takes the play-by-play from a baseball game and turns it into an article in English. Here are various articles about that (written by humans):

They have since widened scope to try and write stories in a huge number of subject domains, and the problems they're trying to solve in changing 'data' into 'stories' sounded both challenging and incredibly interesting. I agreed to an interview.

During this interview, in addition to being subjected to a few hours of questioning (including my first ever attempts at writing code on a white board!), I discovered a lively atmosphere full of fun people working hard on really interesting problems. Furthermore their data team are creating something cool that I think I'll be able to contribute to in a very significant way.

Happilly, they liked me too and I'm really excited to say that I'll be joining the team at Narrative Science on October 30th

## The height and weight of every active football player

While brushing up my Python skills the other week, I decided to figure out web scraping and to do a little data analysis project. Specifically, I decided to look at the height and weight of every active NFL player, with the data that is found in tabular form on nfl.com. For anybody who wants to play with this data for themselves, here is a link to an iPython notebook with the source showing how I made the dataset. This script also requires the file teams.txt, which is up on github via the link. Then, given this database, it is easy to extract the heights and weights of each player and see how they all stack up against one another. A lot of the questions about the graphs were asking about the outlying points, so I'll quickly summarize them here:

• The shortest player is Trindon Holliday, who is 5 feet 5 inches tall
• The 6'8" quarterback is Brock Osweiler

You can click here for a metric version of the chart, translated by the Polish American Football team Angels Toruń! The average weight and height of every position is as follows:

Safeties 6’0″ 207.6
Linebackers 6’2″ 246.3
Defensive Tackles 6’3″ 309.8
Defensive Ends 6’4″ 283.1
Cornerbacks 5’11″ 193.4
Centers 6’3″ 306.2
Tight Ends 6’4″ 254.7
Running Backs 5’11″ 215.3
Guards 6’4″ 314.5
Quarterbacks 6’3″ 223.8
Offensive Tackles 6’5″ 313.5

The same plot for punters and kickers is here. I moved punters and kickers into a separate graph just because the offense and defense are already sufficiently crowded that it's not entirely easy to read.

After posting these visualizations to Reddit they got picked up by the websites of a whole lot of places:

When I add up the number of people who have seen and/or shared these graphs, it's obvious that they have been viewed well over a quarter of a million times (125,000 from a link on imgur.com alone). This is more people than have ever seen the astrophysics that I've worked on for the last decade. I'm not sure how to feel about that.

## Python new-style string formatting reference

I am forever forgetting how to format strings in Python and have probably Googled questions related to Python's format specification mini-language over a thousand times now. For my own reference I put together a quick-reference for the commands I type most frequently. Maybe it will be useful to other chronic Googlers of the same thing, so feel free to use this for anything that you like.

The new-style Python string formatting quick reference is also available as a PDF Ahere

## Crowd Sourcing NFL Predictions #5: Week 3 & 4 Results

Busy couple of weeks, so I got a bit behind with producing these graphs. I'm still very busy so I'll present the graphs here and defer any detailed discussion until later in the season.

### Week 3

This week the games are

• Chargers @ Titans
• Lions @ Redskins
• Giants @ Panthers
• Falcons @ Dolphins
• Bills @ Jets
• Bears @ Steelers

The predictions from the Guardian's Pick Six competition look like this:

Not too bad, 4/6 predicted correctly. And here is how the crowd prediction stacks up against the experts of ESPN, CBS and Yahoo:

Again, the average crowd predicion seems to be holding its own against the experts quite admirably, finishing well above the average expert score this week.

### Week 4

Here are this week's games

• Steelers @ Vikings
• Cardinals @ Buccaneers
• Bears @ Lions
• Seahawks @ Texans
• Cowboys @ Chargers
• Patriots @ Falcons

Along with the Guardian Pick Six aggregate predictions

The Vikings Steelers game was predicted incredibly close, 106 votes to 105 in favour of the Steelers, although you can't see that because the bar is too small. The crowd only called 2/6 of the games correctly this month, and the experts?

A much wider spread, with many of them outperforming the crowd on this week. Weird that the distribution seems so bimodal this week, it would be interesting to check if the experts who predicted 1 correct and 3 correct are clustering around picking AFC and NFC teams or something similar. I have no idea if there was something significant about the games chosen that made this occur, or if it's all in the noise.

### How Good is the Crowd?

I showed this graph with the first few points a couple of weeks ago but have now extended it to include the data through week 4.

There definitely appears to be a positive correlation between the strength with which the crowd leans one way and the score differential, but I still don't want to try and fit anything to it until there is more data. More pressing than the correlation, though, is matplotlib's horrible default color set. I promise to do something about that before presenting further results.

## Crowd Sourcing NFL Predictions #4: Experts vs. The Crowd

I decided to see how, on a person by person basis, the predictions of the experts compared to those made by the general public. On the Guardian Pick Six page, the organizer Paolo Bandini with help from Hamza Mohamed and Andrew Deering posted a comment showing exactly how many games each person predicted correctly. It was easy to turn this comment into a histogram and compare it to the results given by the experts.

Here, the blue bars show the results from the competition entries and the red shows the predictions of the experts. The difference is really striking! We can quantify this by asking what is the probability that the expert scores come from the same parent distribution as the crowd scores with the $\chi^2$ statistic. The answer is that the probability of them coming from the same distribution is vanishingly small.

At least for this first week of results, it appears that the experts are doing worse than the average person. My guess is that because they're being put in public, the experts have a desire to make wild and unlikely predictions so they can grab at the glory if they should somehow come true, although I would be interested to hear anybody else's explanations.

## Crowd Sourcing NFL Predictions #3: Week 2 Results

This is one in my series of posts about comparing NFL predictions made by experts to those sourced from the crowd. In the last post I showed results from six games in week one of the NFL season. In this post, I extend these results to include the second weeks of games. The games that were chosen for the Pick Six competition on the Guardian this week were:

• Miami Dolphins @ Indianapolis Colts
• Carolina Panthers @ Buffalo Bills
• Dallas Cowboys @ Kansas City Chiefs
• Detroit Lions @ Arizona Cardinals
• Denver Broncos @ New York Giants
• San Fransisco 49ers @ Seattle Seahawks

To see how the data was collected, see this post. Let's look at what we find. First, how did the crowdsourced predictions fare? I show here the percentage of predictions that go to each team. Green colored bars show correct predictions, red ones show errors. The crowds didn't do too well

There were only two games for which over 50% of the predictions were correct. A complete turnaround from last week, where the 'average guess' for 5/6 of the games was correct. We don't yet have enough statistics to do anything really interesting, but I'm looking forward to seeing where this goes. Let us now compare the crowd sourced predictions to those made by a sample of 23 experts who are posting predictions publically on ESPN, CBS and Yahoo.

The experts also has a rough week, with the majority of them calling only two out of the six games correctly, although the average expert called 3.08 games correctly.

So, how does the crowd sourced prediction compare to the actual results of the game. The following scatterplot shows the actual score in each game as a function of the percentage of predictions for an away win in each game

The point color shows what week the game is from. A correlation is beginning to appear, but at this point there are only 12 points on the graph and the correlation is largely driven by the red point to the top right, which is the Denver Broncos 41-23 victory over the New York Giants in Week 2.

As we only have two weeks data (12 games) so far, I'm not going to start trying to show how accurate each of the sources of data have turned out to be. That'll have to wait a couple of weeks. For now, I'm happy that this is looking to be an interesting dataset.

## The Naive Bayes Classifier

I've needed to derive and use the naive Bayes classifier a few times over the past years, so for my own reference I'm writing down here all of the steps in the derivation of the classifier along with a couple of example use cases so that I can refresh my memory as needed.

### The Problem

The problem we are solving is as follows: Imagine a training dataset described by $n$ features, $F_1, F_2, ... F_n$ . This training set is classified by hand into one of the classes $C$ . The classifier needs to, given a new set of $n$ features, predict which of the classes $C$ the new data corresponds to. That is, we want to calculate $P(C|F_1\cap F_2\cap ... \cap F_n)$ .

### The Naive Bayesian Probabalistic Model

In the case where there are millions of features, or many values that each feature can take, looking up $P(C|F_1\cap F_2\cap ... \cap F_n)$ in probability tables is not feasible, so we turn to Bayes rule:

Only the numerator of the fraction is interesting for classification, the denominator is just a constant number that depends on the particular set of features chosen, so:

We can make the numerator that remains more tractable by noting that we can break down a bunch of joint probabilities through repeated application of the definition of conditional probability. Concretely, for the case of four variables, $A_1$ through $A_4$ :

More generally, for the case of $n$ variables:

We use this to write for our problem:

Now we need to make the assumption that each of the features $F_i$ is conditionally independent of all of the other features. That is, $P(F_i|C\cap F_j)=P(F_i|C)$ and $P(F_i|C\cap F_j\cap F_l)=P(F_i|C)$ and so on. This removes all of the conditions above and we can simplify the previous equation to

### The Binary Classification Problem

So far we have derived the naive Bayes probability model, we need to add a decision rule to make a classifier. One very common choice is just to pick the most probable class (the maximum a posteriori or MAP decision rule). In general, we do not expect the features to be independent, but the naive Bayesian model is still of a surprising amount of use. In general we don't expect the values for the probabilities of class membership to be accurate if the features are dependent, but as long as they're in the correct order the MAP decision rule will still return the correct.

Imagine now a situation where, for a given set of $n$ features $F_1, F_2, F_3,...,F_n$ we wish to decide whether a particular entry belongs to a class ( $C$ ) or does not ( $\bar{C}$ ). We begin by defining the likelihood ratio:

Here, $P(C)$ and $P(\bar{C})$ represent our knowledge of the prior probabilities. This prior probability distribution might be based on our knowledge of frequencies in the larger population, or on frequency in the training set. It is also important to note that the product is carried out over all features, including those that are negative in the sample we're considering.

If the likelihood ratio is greater than 1, the naive Bayesian classifier predicts that a given entry belongs to class $C$ . In practice, there are frequently many thousands of features and the probabilities of any given feature being present may be very small, so we work in terms of the log likelihood

### Worked Example: Sentiment Analysis

Imagine that you have the following set of messages about a mobile phone

We wish to be able to automatically classify new messages that we receive either as having a positive or negative sentiment. From the overall population, $P(+)=5000/10000=0.5$ and $P(-)=5000/10000=0.5$ . We can now calculate the conditional probabilities for each word being in a positive and negative message.

Note here that we have employed 'smoothing' or 'pseudocounting'. Where a term doesn't appear in a given class (e.g. the word hate appears in no positive messages) we don't want to return a zero as it sends the likelihood calculation haywire. Instead we make it the smallest possible number. Now consider trying to automatically classify the following message

The prefactors outside the fromt of the likelihood ratio cancel out ( $P(+)=P(-)$ ), and multiplying out the remaining products gives

I don't know which way this works out, but it highlights one big weakness with naive Bayes. That is that we have assumed the message is a bag of words, where their relative positions are unimportant. For example, the presence of no in front of pleasure makes this a negative statement that naive Bayes interprets incorrectly.

## Crowd Sourcing NFL Predictions #2: Week 1 Results

### Introduction

This is the second in my series of posts on comparing the accuracy of crowdsourced NFL game predictions from the Guardian's Pick Six competition against the predictions made by experts. This being week 1, I don't have too much data to work with, but here are a couple of graphs.

First, how did the crowdsourced predictions fare? In short, well. Of the six games, the masses selected the winner on 5/6 occasions:

Here, I show for each of the six games, the percentage of the votes that went to each team (the zero in the middle means that the vote was split exactly 50/50). Team names in bold represent the winners, and the bars in green show cases where the aggregate prediction went the same way as the result. The predictions match the actual results with the exception of the one that people felt most strongly about, the Eagles surprise win over the Redskins.

Now, what about the experts? Again, this being the first week we don't have much data, but the following graph shows how each of the 23 experts performed:

The horizontal, red line in this graph shows the crowdsourced result. Taking an average across the experts reveals an average score of 2.95, although I have not yet aggregated together the mean expert prediction to compare that to the crowdsourced result. It is clear so far, at least for this tiny dataset, the crowdsourced results are as good as the very best experts and better than 90% of them.

Oh, and the expert who only got one of six correct was Josh Katzowitz, who apparently makes a living with bold predictions including today's article about how the (2-14 last year) Chiefs are going to surprise the (8-8 last year) Cowboys in week 2

## Crowd Sourcing NFL Predictions #1: Data Collection

### Introduction

"The Wisdom of Crowds" is a concept that, since its publication in 2004, has become part of the popular imagination. In short, the wisdom of crowds is the idea that if you aggregate together a large crowd of people, the average decisions they reach in response to difficult questions are likely to be better than those given by any one individual. This comes about because each member of the crowd brings their own viewpoint and set of information to the question, and the crowd as a whole has a more complete and diverse set of information about the problem at hand than does any of its members. The quintessential example of this is the tale of a crowd of people being asked to guess how many jellybeans are in a large jar, or guess the weight of a cow. For large crowds, the average guess is usually closer to the correct answer than any one person's estimate. Indeed, I even remember doing this experiment at school as a child in school assembly and seeing it work first hand.

The start of the football season got me thinking about this, especially given the huge number of football predictions that are available on the internet. Of course, not all crowds are 'wise'. Investors in a stock market bubble, for example, all make the same errors. But it did strike me that crowdsourced NFL predictions probably do meet all of the criteria necessary for a crowd to make very good judgments (cribbed from the Wikipedia page on The Wisdom of Crowds:

• Diversity of opinion (i.e. Each person should have private information even if it's just an eccentric interpretation of the known facts.)
• Independence (i.e. people's opinions are not fixed by those of their neighbors)
• Decentralization (i.e. people can draw on local knowledge)
• Aggregation (i.e. Some mechanism exists for turning private judgments into a collective decision.)

This certainly sounds a lot like football predictions on the Internet. People from the world over have different sets of information and prejudices about their teams, they're sitting alone in front of a computer, and we can easily pull together a bunch of predictions to provide an aggregate prediction for the outcome of games.

I decided to run a little experiment. For the 2013 NFL season, I will compare the accuracy of a crowd of people predicting NFL game results (The Guardian's Pick Six competition) to the set of experts on ESPN, Fox and the like. I describe below where each dataset comes from, then in future posts I'll show the results.

### Crowd Sourced Data

This data comes from The Guardian's NFL Pick Six competition. This is a bit of fun on the website of the Guardian newspaper where six games are selected each week, and the readers are invited to make their predictions for the winner of each one. I like using this dataset because it immediately cuts out the more interesting matchups each week. The organizer of Pick Six selects the games where the outcome is less than obvious.

Getting the data from the website is a non-trivial process. People submit their predictions as comments on the article so we need to do a significant amount of cleaning to retreive the data. Here are a few examples of prediction posts from week 1:

This is not a homogeneous dataset. Some people list team names, others list cities, some leave comments or jokes, and we need to worry about typos. This is a job that a human can deal with very easily. Indeed, I tallied up the first 100 (of 700 total) predictions by hand, but it's a time consuming job that I do not have time to do every week. Here are the results for the first 100 or so predictions:

We need, therefore, to find an automated way of getting the predictions from the posts. I extracted an approximate set of predictions automatically with the following steps:

• Take a copy of the complete text of all of the comments from the webpage
• Extract and count up the number of lines that match the following criteria
• Line begins with the team name (case insensitive), AND
• Line does not include the character '@'

This is all encapsulated in this simple command line

cat week1.dat | grep -i ^teamname | grep -v '@' | wc -l

The reason for removing the @ symbol is that a number of people write predictions in the form Bengals@Bears, and use boldface text to highlight the winner. Without removing the lines with @ symbols, this simple method assigns all predictions formatted in this way to the team listed first (the away team). This is obviously simplified method that makes some errors. For example, it misses anybody who uses the city as a proxy for teamname, or makes typos, or writes variations on teamnames (e.g. a number of people purposely use 'Washington' or a euphemism in place of 'Redskins' because of the controversy surrounding this name). Nevertheless, over 85% of posts are of a form (team names, one per line) that can easily be handled by this script.

To check how accurate the results from this simple data scraper are, I compared its results with the ones that I obtained doing everything by hand. In the following graph, I show the percentage of people who predict a victory for the away team for the exact (i.e. where I read every post and tally the correct result) and the approximate (i.e. I ran the command line above to extract the number of predictions) methods

The two always agree to within a percentage point or two, So we now have an automated method of extracting how many Guardian Pick Six predictions favor each team. Throughout the rest of this series I will be using this method on every prediction posted before the kickoff of the relevant games.

### Data from "Experts"

A number of news outlets publish NFL predictions. I'll be taking the picks from the following experts: First, ESPN provides picks from 13 NFL experts at this website in tabular form:

CBS provides picks from another 8 experts at this website in tabular form:

Finally, Yahoo provides picks from 2 experts at this website in tabular form:

This gives us a total sample of 23 NFL experts, and I'll also aggregate the experts together into a single 'expert prediction'. Results to be posted soon after the Monday night football finishes tonight.

## Analysis of 250,000 geotagged tweets from Chicago

The City of Chicago provides a rather comprehensive data portal, full of datasets related to the city. I wanted to take a look at these as an excuse to learn some R, a language I don't have too much experience with as yet. I decided that on top of the official city dataset, I would add one of my own and as such wrote (in Python) a robot that watches the Twitter public stream for geotagged tweets from the Chicago area. Whenever it finds one it stores the following in an SQL database:

• The tweet text, username and date
• The GPS coordinates from which the tweet was sent
• Which of Chicago's 77 community areas the tweet was sent from
• Whether or not the tweet has a positive sentiment

The neighborhood that a tweet comes from is calculated by comparing the GPS coordinates to the shapefiles that define the edges of Chicago's neighborhoods, as supplied by the City of Chicago. The final quantity, sentiment, is calculated in a simple way. I just compare every word in the tweet against a wordlist that contains for 2,500 common words an estimate of how positive or negative each one is (full disclosure: I got this list as part of an assignment while taking Coursera's Data Science course). A few examples of words in here:

abhors -3
abilities 2
attraction 2
awesome 4
awful -3
awkward -2

The happiness of each tweet is then simply the happiness of each of its constituent words. This is obviously not a well calibrated method, but it should allow us to compare different populations of tweets.

I ran this bot for a full week in May of 2013 and ended up collecting just over a quarter of a million tweets to analyse. Here are a few examples of happy tweets (anonymized):

Had fun great time with my boyfriend... Time for bed.. I hope every had a great night because I did.. Good Night..

Good morning lake. Good morning runners. Good morning coffee shop. Good morning hipsters. Good morning Chicago.

@anonymous I love you and admire you. Thanks for making Louis happy. Have a lovely day. #HappyBirthday

I don't really feel comfortable sharing people's most upset tweets, even if they are anonymized, as there is some surprisingly personal material in there (along with a LOT of swearing). We can then put all of this together to see which of Chicago's neighborhoods have the happiest tweets. I don't think that the results will be much of a surprise to anybody who lives here:

The well-off North of the city seems to make systematically happier tweets than the neighborhoods to the south. Additionally, you can't see this from the map, but 100% of the tweets sent from on the water were happy, which I guess is a good argument for getting a boat. Anyway, this map got me interested in asking how the economic status of a neighborhood affected its conversation, and I show this here:

Here, each point is a different Chicago neighborhood. The size of each point is proportional to the number of tweets that were sent, the color of each point is proportional to the number of violent crimes per person (purple is low, yellow is high), and the axes show the average happiness of the tweets as a function of the local unemployment rate.

The correlation between average mood and unemployment rate is strikingly strong, with the exception of one neighborhood to the top-right of the graph. This is Garfield Park, and on further investigation I found that it was one single person who tweeted "Good Night, @anonymous, My Love" to over 200 people, accounting for a good fraction of the neighborhood's total tweets.

This is a pretty cool dataset, and this is only a very preliminary analysis, but with many thousands of words per neighborhood I plan to look at how word usage clusters around neighborhoods of different types at som epoint in the future.


text:
# What is the number of roots of this equation?

We are given ordered sets of positive integers A and B with n elements each and an integer C, how can we find the number of integer roots of the following equation:
$\sum_{i=1}^n floor(\frac{A_i*x}{B_i}) = C$
I have been stuck on this problem for such a long time. I have tried solving using all the properties of floor functions I know till now, it just doesn't fit into any.
An example will be when A = (1,5) and B = (2,4) and C = 3, the number of roots is 1, the root is 5. Any help will be very helpful to me, thanks.

• Welcome to Math.SE Kindly share your thoughts on the question you put up so as to maximize the number of possible contributors. – The Dead Legend Apr 28 '17 at 8:39
• Done :) Thanks for trying to help me out. I am new here. – Shreya Halder Apr 28 '17 at 9:11

If $x=x_0$ is a solution then so is any x in a small enough interval to the right or the left of $x_0$, so an equation of this form will not have discrete roots - its solution set will consist of one or more intervals (unless it is empty).
If all the rationals $A_i/B_i$ are positive then the sum of the floor functions is a sum of monotonically increasing functions so it is itself a monotonically increasing function of x. Its solution set is therefore either empty or it consists of a single interval. An example with an empty solution set is A=(1,1), B=(1,1), C=1 i.e. $$2 floor(x) = 1$$
Similar argument holds if all $A_i/B_i$ are negative.
If $A_i/B_i$ are both positive and negative then the solution set can consist of an infinite number of intervals. For example, consider A=(-1,-1,2), B=(1,1,1), C=0 i.e.
$$2floor(-x) +floor(2x)=0$$


text:
Quantitative comparisons of macrophotography with and without Lieberkühn reflectors

In order to quantitatively examine the effect of the the 3D printed Lieberkühn reflectors I described previously, I came up with two image quality metrics relevant to their use, both measured on the “dark side” of the image subject. The metrics we will look at today are average intensity and, as a measure of contrast, the standard deviation of pixel values on a line trace.

I’ll be using the Megachile photo from the previous post for these analyses.

The first line trace begins at right eye and extends back behind the wing:

If we plot these values together, we see that the photo taken with the Lieberkühn (values in black) is brighter and brings out a lot more detail, while the photo taken without is relatively flat and dark. We see similar results for second and third traces, taken across the tegula and along the wing.

If we compare the average values:

octave3.2:25> mean(RE523(:,2)) %with Lieberkühn, right eye trace
ans = 102.00
octave3.2:26> mean(RE524(:,2))%w/o Lieberkühn, right eye trace
ans = 54.093
octave3.2:28> mean(AT523(:,2))%with Lieberkühn, trace across tegula
ans = 81.553
octave3.2:27> mean(AT524(:,2))%w/o Lieberkühn, trace across tegula
ans = 51.288
octave3.2:29> mean(W523(:,2))%with Lieberkühn, across the wing
ans = 103.85
octave3.2:30> mean(W524(:,2))%w/o Lieberkühn, across the wing
ans = 53.045

We see that taken together, the averages of the plots from the photo taken with the Lieberkühn are about 80% brighter than those without.

mean(Lieberkühn)/mean(no Lieberkühn) = 1.8142

octave3.2:56> std(RE523(:,2)))%with Lieberkühn, right eye trace
ans = 20.316
octave3.2:55> std(RE524(:,2))%w/o Lieberkühn, right eye trace
ans = 7.3926

octave3.2:54> std(AT523(:,2))%with Lieberkühn, trace across tegula
ans = 17.737
octave3.2:53> std(AT524(:,2))%w/o Lieberkühn, trace across tegula
ans = 13.227

octave3.2:52> std(W523(:,2))%with Lieberkühn, across the wing
ans = 12.746
octave3.2:51> std(W524(:,2))%w/o Lieberkühn, across the wing
ans = 8.2902

Using standard deviation as a metric for image detail, we get an increase of about 75% in standard deviation over the dark photo by using the reflector.

octave3.2:30> (20.316+17.737+12.746)/(7.3926+13.227+8.2902)
ans = 1.7572

The averages, standard deviation etc. may seem a bit redundant at this point; you don’t need to plot a pixel-value profile to see that the image with the reflector is much brighter and more detailed than the photo taken without.

3D Printable Lens Hood Design

A lens hood is a shade that blocks out-of-frame light from reflecting off of the internals within the lens. This minimizes lens flares, so you can add them later in post. Just kidding.

Another form of lens flare is less obvious (and I don’t think J.J. Abrams uses it). It manifests as a haze across the majority of the frame making the image appear washed-out, and it never looks good. Unlike deliberate lens flares, it’s not obvious in the image itself where it comes from and doesn’t look dramatic.

To get the most effect from a lens hood, it needs to block out as much unwanted light as possible without actually showing up in the frame. This means that for any given lens at a certain focal length and field of view there will be a best angle for your lens hood.

The wikipedia article for angle of view gives an equation depending on the focal length and sensor size.

$2cdot tan^{-1}(frac{d}{2f})$

The variable $d$ is the dimension of interest. For a lens hood with a simple circular cross section throughout the longest dimension should be used, e.g. the diagonal length of a typical rectilinear sensor. The doubling factor can be omitted if you want to work the angle in relation to the optical axis, rather than the total angle.

The lens hood below is a general purpose lens hood (also 3D printed) for lenses with a 58mm filter thread diameter. It flares out a bit, and the angle is wide enough to use with a ~27mm focal length lens.

The images below show essentially the same 58mm diameter lens hood optimized for 16mm, 35mm, and 50mm, in order from left to right. The length of the hood in each case is 16mm. The shorter the focal length of the lens (and the larger the image sensor) the wider the angle, and the lens hood angle increases accordingly.

So far, I have printed the general purpose lens hood, which errs on the side of wide-angle caution. Once I have the additional test pieces in hand, we’ll give ’em the old Pepsi challenge.


text:
# Prove that there are infinitely many natural numbers that can't be written as $a^2+p$

Generally speaking, if they ask us to prove that there are infinitely many numbers that can't be written in a certain way, how should we try to solve the problem? I've never seen a solution to such a problem yet, so I can't think of any strategy to attack the problem in the title. I think the problem leaves too many possibilities too check! If they tell me that an integer number can't be written in the form $a^2+p$ where $p$ is prime that doesn't give me any clues about how such numbers could be. I mean it doesn't put enough restrictions on the numbers that aren't of that form that I go after studying the properties of such numbers.

There is also another problem of the similar type in my book: Prove that there are infinitely many odd numbers such that they can't be written as the sum of less than 3 prime numbers.

-
The other problem is easy, if an odd number can be written as the sum of two primes, one of the two primes must be even. That restricts the possibilities considerably. –  Daniel Fischer Aug 30 '13 at 21:50

Think of numbers that are perfect squares. So is it possible to write a perfect square of $x$ as:$$x^2=a^2+p$$

-

Try to think of what kinds of numbers might have restrictions on how they can be written, in terms of primes and squares etc. since your problem involves primes and squares. Hint: What if $n = q^2$ where $q$ is an integer?

-

For your first problem, one maybe looks for numbers that are "unlikely" to be of the desired form. Take for example a perfect square $b^2$. In order for $b^2$ to be of shape $a^2+p$, we must have $b^2-a^2=p$. But because of the natural factorization $(b-a)(b+a)$, that forces $a=b-1$. So $b+a=2b-1$ must be prime.

So take any number odd $2b-1$ which is not prime. There are infinitely many of these, for example $3^2,3^3,3^4,\dots$.

Then $b^2$ cannot be expressed as $a^2+p$.

-

If $b^2=a^2+p$ then $b=a+1$ and $2a+1=p$. Now consider the numbers of the form $(3k+2)^2$.

-
Or, $(b-a)(b+a)=p^2$ and force $p$ not to be prime :) –  chubakueno Sep 1 '13 at 21:50

The proof is given in the OEIS at A014090, where it's noted that $21679$ is the last known non-square not of the form $a^2+p$. The comments and references at A020495 are also worth a look.

-

For your second problem, "less than three" means "one or two".

So, what odd numbers can't be a prime (i.e., sum of one prime) or the sum of two primes?

Is there a special prime which might help solve the second part (sum of two primes)?

-


text:
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_rand_bb_inc (g05xd)

## Purpose

nag_rand_bb_inc (g05xd) computes scaled increments of sample paths of a free or non-free Wiener process, where the sample paths are constructed by a Brownian bridge algorithm. The initialization function nag_rand_bb_inc_init (g05xc) must be called prior to the first call to nag_rand_bb_inc (g05xd).

## Syntax

[z, b, ifail] = g05xd(npaths, diff, z, c, rcomm, 'rcord', rcord, 'd', d, 'a', a)
[z, b, ifail] = nag_rand_bb_inc(npaths, diff, z, c, rcomm, 'rcord', rcord, 'd', d, 'a', a)

## Description

For details on the Brownian bridge algorithm and the bridge construction order see Section [Brownian Bridge] in the G05 Chapter Introduction and Section [Description] in (g05xc). Recall that the terms Wiener process (or free Wiener process) and Brownian motion are often used interchangeably, while a non-free Wiener process (also known as a Brownian bridge process) refers to a process which is forced to terminate at a given point.
Fix two times t0 < T${t}_{0}, let (ti)1in ${\left({t}_{i}\right)}_{1\le i\le \mathit{n}}$ be any set of time points satisfying t0 < t1 < t2 < < tn < T${t}_{0}<{t}_{1}<{t}_{2}<\cdots <{t}_{\mathit{n}}, and let Xt0${X}_{{t}_{0}}$, (X ti )1in ${\left({X}_{{t}_{i}}\right)}_{1\le i\le \mathit{n}}$, XT${X}_{T}$ denote a d$d$-dimensional Wiener sample path at these time points.
The Brownian bridge increments generator uses the Brownian bridge algorithm to construct sample paths for the (free or non-free) Wiener process X$X$, and then uses this to compute the scaled Wiener increments
 ( Xt1 − Xt0 )/( t1 − t0 ) , ( Xt2 − Xt1 )/( t2 − t1 ) , … , ( Xtn − Xtn − 1 )/( tn − tn − 1 ) , ( XT − Xtn )/( T − tn ) $Xt1 - Xt0 t1 - t0 , Xt2 - Xt1 t2 - t1 ,…, Xtn - Xtn-1 tn - tn-1 , XT - Xtn T - tn$
The example program in Section [Example] shows how these increments can be used to compute a numerical solution to a stochastic differential equation (SDE) driven by a (free or non-free) Wiener process.

## References

Glasserman P (2004) Monte Carlo Methods in Financial Engineering Springer

## Parameters

Note: the following variable is used in the parameter descriptions: n = ntimes$\mathit{n}={\mathbf{ntimes}}$, the length of the array times passed to the initialization function nag_rand_bb_inc_init (g05xc).

### Compulsory Input Parameters

1:     npaths – int64int32nag_int scalar
The number of Wiener sample paths.
Constraint: npaths1${\mathbf{npaths}}\ge 1$.
2:     diff(d) – double array
d, the dimension of the array, must satisfy the constraint d1${\mathbf{d}}\ge 1$.
The difference between the terminal value and starting value of the Wiener process. If a = 0${\mathbf{a}}=0$, diff is ignored.
3:     z(ldz, : $:$) – double array
The first dimension, ldz, of the array z must satisfy
• if rcord = 1${\mathbf{rcord}}=1$, ldzd × (n + 1a)$\mathit{ldz}\ge {\mathbf{d}}×\left(\mathit{n}+1-{\mathbf{a}}\right)$;
• if rcord = 2${\mathbf{rcord}}=2$, ldznpaths$\mathit{ldz}\ge {\mathbf{npaths}}$.
The second dimension of the array must be at least npaths${\mathbf{npaths}}$ if rcord = 1${\mathbf{rcord}}=1$ and at least d × (n + 1a)${\mathbf{d}}×\left(\mathit{n}+1-{\mathbf{a}}\right)$ if rcord = 2${\mathbf{rcord}}=2$
The Normal random numbers used to construct the sample paths.
If quasi-random numbers are used, the d × (n + 1a)${\mathbf{d}}×\left(\mathit{n}+1-{\mathbf{a}}\right)$-dimensional quasi-random points should be stored in successive rows of Z$Z$.
4:     c(ldc, : $:$) – double array
The first dimension of the array c must be at least d${\mathbf{d}}$
The second dimension of the array must be at least d${\mathbf{d}}$
The lower triangular Cholesky factorization C$C$ such that CCT$C{C}^{\mathrm{T}}$ gives the covariance matrix of the Wiener process. Elements of c above the diagonal are not referenced.
5:     rcomm(12 × (ntimes + 1)$12×\left(\mathit{ntimes}+1\right)$) – double array
Communication array as returned by the last call to nag_rand_bb_inc_init (g05xc) or nag_rand_bb_inc (g05xd). This array must not be directly modified.

### Optional Input Parameters

1:     rcord – int64int32nag_int scalar
The order in which Normal random numbers are stored in z and in which the generated values are returned in b.
Default: 1$1$
Constraint: rcord = 1${\mathbf{rcord}}=1$ or 2$2$.
2:     d – int64int32nag_int scalar
Default: The dimension of the array diff and the first dimension of the array c. (An error is raised if these dimensions are not equal.)
The dimension of each Wiener sample path.
Constraint: d1${\mathbf{d}}\ge 1$.
3:     a – int64int32nag_int scalar
If a = 0${\mathbf{a}}=0$, a free Wiener process is created and diff is ignored.
If a = 1${\mathbf{a}}=1$, a non-free Wiener process is created where diff is the difference between the terminal value and the starting value of the process.
Default: 1$1$
Constraint: a = 0${\mathbf{a}}=0$ or 1$1$.

ldz ldc ldb

### Output Parameters

1:     z(ldz, : $:$) – double array
The first dimension, ldz, of the array z will be
• if rcord = 1${\mathbf{rcord}}=1$, ldzd × (n + 1a)$\mathit{ldz}\ge {\mathbf{d}}×\left(\mathit{n}+1-{\mathbf{a}}\right)$;
• if rcord = 2${\mathbf{rcord}}=2$, ldznpaths$\mathit{ldz}\ge {\mathbf{npaths}}$.
The second dimension of the array will be npaths${\mathbf{npaths}}$ if rcord = 1${\mathbf{rcord}}=1$ and at least d × (n + 1a)${\mathbf{d}}×\left(\mathit{n}+1-{\mathbf{a}}\right)$ if rcord = 2${\mathbf{rcord}}=2$
The Normal random numbers premultiplied by c.
2:     b(ldb, : $:$) – double array
The first dimension, ldb, of the array b will be
• if rcord = 1${\mathbf{rcord}}=1$, ldbd × (n + 1)$\mathit{ldb}\ge {\mathbf{d}}×\left(\mathit{n}+1\right)$;
• if rcord = 2${\mathbf{rcord}}=2$, ldbnpaths$\mathit{ldb}\ge {\mathbf{npaths}}$.
The second dimension of the array will be npaths${\mathbf{npaths}}$ if rcord = 1${\mathbf{rcord}}=1$ and at least d × (n + 1)${\mathbf{d}}×\left(\mathit{n}+1\right)$ if rcord = 2${\mathbf{rcord}}=2$
The scaled Wiener increments.
Let Xp,ik${X}_{p,i}^{k}$ denote the k$k$th dimension of the i$i$th point of the p$p$th sample path where 1kd$1\le k\le {\mathbf{d}}$, 1in + 1$1\le i\le \mathit{n}+1$ and 1pnpaths$1\le p\le {\mathbf{npaths}}$. The increment ((Xp,ikXp,i1k))/((titi1))$\frac{\left({X}_{p,i}^{k}-{X}_{p,i-1}^{k}\right)}{\left({t}_{i}-{t}_{i-1}\right)}$ is stored at B(p,k + (i1) × d)$B\left(p,k+\left(i-1\right)×{\mathbf{d}}\right)$.
3:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
On entry, rcomm was not initialized or has been corrupted. On entry, rcomm was not initialized or has been corrupted. On entry, rcomm was not initialized or has been corrupted.
ifail = 2${\mathbf{ifail}}=2$
Constraint: npaths1${\mathbf{npaths}}\ge 1$.
ifail = 3${\mathbf{ifail}}=3$
On entry, the value of rcord is invalid.
ifail = 4${\mathbf{ifail}}=4$
Constraint: d1${\mathbf{d}}\ge 1$.
ifail = 5${\mathbf{ifail}}=5$
Constraint: a = 0 or 1${\mathbf{a}}=0\text{ or }1$.
ifail = 6${\mathbf{ifail}}=6$
Constraint: ldzd × (ntimes + 1a)$\mathit{ldz}\ge {\mathbf{d}}×\left({\mathbf{ntimes}}+1-{\mathbf{a}}\right)$.
Constraint: ldznpaths$\mathit{ldz}\ge {\mathbf{npaths}}$.
ifail = 7${\mathbf{ifail}}=7$
ldc is too small.
ifail = 8${\mathbf{ifail}}=8$
Constraint: ldbd × (ntimes + 1)$\mathit{ldb}\ge {\mathbf{d}}×\left({\mathbf{ntimes}}+1\right)$.
Constraint: ldbnpaths$\mathit{ldb}\ge {\mathbf{npaths}}$.
ifail = 999${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

Not applicable.

None.

## Example

```function nag_rand_bb_inc_example
% We wish to solve the stochastic differential equation (SDE)
%   dSt = r*St*dt + sigma*St*dXt
% where X is a one dimensional Wiener process.
% This means we have
%     a = 0
%     d = 1
%     c = 1
% We now set the other parameters of the SDE and the Euler-Maruyama scheme

% Initial value of the process
s0 = 1;
r = 0.05;
sigma = 0.12;

% Number of paths to simulate
npaths = int64(3);

% The time interval [t0,T] on which to solve the SDE
t0   = 0;
tend = 1;

% The time steps in the discretization of [t0,T]
ntimesteps = 20;
t = t0 + (1:ntimesteps)*(tend-t0)/(ntimesteps+1);

% Make the bridge construction order
bgord = int64(3);
move = zeros(0, 1, 'int64');
[times, ifail] = nag_rand_bb_make_bridge_order(t0, tend, t, move, 'bgord', bgord);

% Generate the Z values
d = 1;
a = int64(0);
[z] = get_z(npaths, d, a, ntimesteps);

% Initialize the generator
[rcomm, ifail] = nag_rand_bb_inc_init(t0, tend, times);

% Get the scaled increments of the Wiener process
dif = zeros(d, 1);
c   = ones(d, 1);
[z, b, ifail] = nag_rand_bb_inc(npaths, dif, z, c, rcomm, 'a', a);

% Do the Euler-Maruyama time stepping
st = zeros(ntimesteps+1, npaths);

% Do first time step
st(1,:) = s0 + (t(1)-t0)*(r*s0+sigma*s0*b(1, :));
for i=2:ntimesteps
st(i,:) = st(i-1,:) + (t(i)-t(i-1))*(r*st(i-1,:)+sigma*st(i-1,:).*b(i,:));
end
% Do last time step
st(ntimesteps+1,:) = st(ntimesteps,:) + (tend-t(ntimesteps))*...
(r*st(ntimesteps,:)+sigma*st(ntimesteps,:).*b(ntimesteps+1,:));

% Compute the analytic solution:
%    ST = S0*exp( (r-sigma^2/2)T + sigma WT )
analytic = s0*exp((r-0.5*sigma*sigma)*tend+sigma*sqrt(tend-t0)*z(1,:));

% Display the results
fprintf('\nEuler-Maruyama solution for Geometric Brownian motion\n');
fprintf('        Path 1     Path 2     Path 3\n');
for i = 1:ntimesteps+1
fprintf('%2d  %10.4f %10.4f %10.4f\n', i, st(i, :));
end
fprintf('\nAnalytic solution at final time step\n');
fprintf('        Path 1     Path 2     Path 3\n');
fprintf('    %10.4f %10.4f %10.4f\n', analytic);

function [z] = get_z(npaths, d, a, ntimes)
idim = d*(ntimes+1-a);

% We now need to generate the input pseudorandom points

% First initialize the base pseudorandom number generator
state = initialize_prng(int64(6), int64(0), [int64(1023401)]);

% Scrambled quasi-random sequences preserve the good discrepancy
% properties of quasi-random sequences while counteracting the bias
% some applications experience when using quasi-random sequences.
% Initialize the scrambled quasi-random generator.
[iref, state] = initialize_scrambled_qrng(int64(1), int64(2), idim, state);

% Generate the quasi-random points from N(0,1)
xmean = zeros(idim, 1);
std   = ones(idim, 1);
[z, iref, ifail] = nag_rand_quasi_normal(xmean, std, npaths, iref);
z = z';
function [state] = initialize_prng(genid, subid, seed)
% Initialize the generator to a repeatable sequence
[state, ifail] = nag_rand_init_repeat(genid, subid, seed);

function [iref, state] = initialize_scrambled_qrng(genid,stype,idim,state)
iskip = int64(0);
nsdigits = int64(32);
[iref, state, ifail] = nag_rand_quasi_init_scrambled(genid, stype, int64(idim), iskip, nsdigits, state);
```
```

Euler-Maruyama solution for Geometric Brownian motion
Path 1     Path 2     Path 3
1      0.9668     1.0367     0.9992
2      0.9717     1.0254     1.0077
3      0.9954     1.0333     1.0098
4      0.9486     1.0226     0.9911
5      0.9270     1.0113     1.0630
6      0.8997     1.0127     1.0164
7      0.8955     1.0138     1.0771
8      0.8953     0.9953     1.0691
9      0.8489     1.0462     1.0484
10      0.8449     1.0592     1.0429
11      0.8158     1.0233     1.0625
12      0.7997     1.0384     1.0729
13      0.8025     1.0138     1.0725
14      0.8187     1.0499     1.0554
15      0.8270     1.0459     1.0529
16      0.7914     1.0294     1.0783
17      0.8076     1.0224     1.0943
18      0.8208     1.0359     1.0773
19      0.8190     1.0326     1.0857
20      0.8217     1.0326     1.1095
21      0.8084     0.9695     1.1389

Analytic solution at final time step
Path 1     Path 2     Path 3
0.8079     0.9685     1.1389

```
```function g05xd_example
% We wish to solve the stochastic differential equation (SDE)
%   dSt = r*St*dt + sigma*St*dXt
% where X is a one dimensional Wiener process.
% This means we have
%     a = 0
%     d = 1
%     c = 1
% We now set the other parameters of the SDE and the Euler-Maruyama scheme

% Initial value of the process
s0 = 1;
r = 0.05;
sigma = 0.12;

% Number of paths to simulate
npaths = int64(3);

% The time interval [t0,T] on which to solve the SDE
t0   = 0;
tend = 1;

% The time steps in the discretization of [t0,T]
ntimesteps = 20;
t = t0 + (1:ntimesteps)*(tend-t0)/(ntimesteps+1);

% Make the bridge construction order
bgord = int64(3);
move = zeros(0, 1, 'int64');
[times, ifail] = g05xe(t0, tend, t, move, 'bgord', bgord);

% Generate the Z values
d = 1;
a = int64(0);
[z] = get_z(npaths, d, a, ntimesteps);

% Initialize the generator
[rcomm, ifail] = g05xc(t0, tend, times);

% Get the scaled increments of the Wiener process
dif = zeros(d, 1);
c   = ones(d, 1);
[z, b, ifail] = g05xd(npaths, dif, z, c, rcomm, 'a', a);

% Do the Euler-Maruyama time stepping
st = zeros(ntimesteps+1, npaths);

% Do first time step
st(1,:) = s0 + (t(1)-t0)*(r*s0+sigma*s0*b(1, :));
for i=2:ntimesteps
st(i,:) = st(i-1,:) + (t(i)-t(i-1))*(r*st(i-1,:)+sigma*st(i-1,:).*b(i,:));
end
% Do last time step
st(ntimesteps+1,:) = st(ntimesteps,:) + (tend-t(ntimesteps))*...
(r*st(ntimesteps,:)+sigma*st(ntimesteps,:).*b(ntimesteps+1,:));

% Compute the analytic solution:
%    ST = S0*exp( (r-sigma^2/2)T + sigma WT )
analytic = s0*exp((r-0.5*sigma*sigma)*tend+sigma*sqrt(tend-t0)*z(1,:));

% Display the results
fprintf('\nEuler-Maruyama solution for Geometric Brownian motion\n');
fprintf('        Path 1     Path 2     Path 3\n');
for i = 1:ntimesteps+1
fprintf('%2d  %10.4f %10.4f %10.4f\n', i, st(i, :));
end
fprintf('\nAnalytic solution at final time step\n');
fprintf('        Path 1     Path 2     Path 3\n');
fprintf('    %10.4f %10.4f %10.4f\n', analytic);

function [z] = get_z(npaths, d, a, ntimes)
idim = d*(ntimes+1-a);

% We now need to generate the input pseudorandom points

% First initialize the base pseudorandom number generator
state = initialize_prng(int64(6), int64(0), [int64(1023401)]);

% Scrambled quasi-random sequences preserve the good discrepancy
% properties of quasi-random sequences while counteracting the bias
% some applications experience when using quasi-random sequences.
% Initialize the scrambled quasi-random generator.
[iref, state] = initialize_scrambled_qrng(int64(1), int64(2), idim, state);

% Generate the quasi-random points from N(0,1)
xmean = zeros(idim, 1);
std   = ones(idim, 1);
[z, iref, ifail] = g05yj(xmean, std, npaths, iref);
z = z';
function [state] = initialize_prng(genid, subid, seed)
% Initialize the generator to a repeatable sequence
[state, ifail] = g05kf(genid, subid, seed);

function [iref, state] = initialize_scrambled_qrng(genid,stype,idim,state)
iskip = int64(0);
nsdigits = int64(32);
[iref, state, ifail] = g05yn(genid, stype, int64(idim), iskip, nsdigits, state);
```
```

Euler-Maruyama solution for Geometric Brownian motion
Path 1     Path 2     Path 3
1      0.9668     1.0367     0.9992
2      0.9717     1.0254     1.0077
3      0.9954     1.0333     1.0098
4      0.9486     1.0226     0.9911
5      0.9270     1.0113     1.0630
6      0.8997     1.0127     1.0164
7      0.8955     1.0138     1.0771
8      0.8953     0.9953     1.0691
9      0.8489     1.0462     1.0484
10      0.8449     1.0592     1.0429
11      0.8158     1.0233     1.0625
12      0.7997     1.0384     1.0729
13      0.8025     1.0138     1.0725
14      0.8187     1.0499     1.0554
15      0.8270     1.0459     1.0529
16      0.7914     1.0294     1.0783
17      0.8076     1.0224     1.0943
18      0.8208     1.0359     1.0773
19      0.8190     1.0326     1.0857
20      0.8217     1.0326     1.1095
21      0.8084     0.9695     1.1389

Analytic solution at final time step
Path 1     Path 2     Path 3
0.8079     0.9685     1.1389

```


text:
# Permutation tests for paired data with several variables

I am trying to answer my own previous question, with a permutation method.

I have a reference method $$M_{ref}$$ that I would like to compare with other methods $$M_1$$, $$M_2$$, $$M_3$$ when taking measurements (denoted scores from now on) from a data sample of size $$N (N > 30)$$. However, when comparing $$M_{ref}$$ with a method $$M_i$$, i do not have one kind of continuous score $$S_j$$, but several, say $$S_1$$, $$S_2$$, $$S_3$$, and $$S_4$$. So for each method, I have $$N$$ scores of type $$S1$$, $$N$$ scores of type $$S2$$, etc.

My question (still) is: how to obtain a non-parametric test p-value for each comparison $$C_{i,j} = _{score=C_j}$$. Please note I am not interested in any comparison between methods other than the reference method (e.g. $$M_{2}$$ Vs $$M_3$$ is of no interest), and comparisons between different kind of scores would not make sense.

So now i am trying to think of a way to test all the scores at the same time with permutation tests for a given method comparison, e.g. $$M_{ref}$$ Vs $$M1$$, and then simply use Bonferonni correction by multiplying by the number of times I compare methods this way (so by 4, because there are 4 methods comparison to do: $$M_{ref}$$ Vs $$M1$$, $$M_{ref}$$ Vs $$M2$$, etc.).

Now I am not sure of the proper way to test all the scores at the same time with permutation tests. Would the following (inspired from here) work?

    # Number of random reps
nreps <- 10000
# All these operations are realized across score columns,
# but randomly assigning sign the same way for all columns
# [is that the main trick?]
# Now we will randomly assign the sign of the difference
Mref.vs.M1.S1.to.S4 <- Mref.S1.to.S4 - M1.S1.to.S4
sample.means <- mean(Mref.vs.M1.S1.to.S4)
cat("Mean differences for S1, S2, S3 and S4 are", sample.means, "\n")
unsigned.diff <- abs(Mref.vs.M1.S1.to.S4)
n <- length(unsigned.diff)
for (i in 1:nreps) {
# Draw N obs randomly from 1, -1
signs <- sample(c(1, -1), size = n, replace = TRUE)
signed.diff <- unsigned.diff*signs
# Calculate the mean with that arrangement of signs
# for S1, S2, S3 and S4
means.random[i] <- mean(signed.diff)
}

# The question is what percentage of those differences
# are greater than or equal to plus or minus sample.means
# of S1, S2, S3 and S4
percent.above <- length(means.random[means.random >= abs(sample.measn)])/nreps
percent.below <- length(means.random[means.random <= abs(sample.means)*(-1)])/nreps
percent.extreme <- percent.above + percent.below

cat("The randomization test gives a probability under the null")
cat("for S1, S2, S3 and S4 [4 p-values] of", percent.extreme, "\n\n")

Does this work or is it just equivalent to doing resampling multiple times for each score type, and the resulting p-values still need to be Bonferonni-corrected by multiplying by the number of score types (4 in this case) ?


text:
apollo_classAlloc {apollo} R Documentation

## Calculates class allocation probabilities for a Latent Class model

### Description

Calculates class allocation probabilities for a Latent Class model using a Multinomial Logit model and can also perform other operations based on the value of the functionality argument.

### Usage

apollo_classAlloc(classAlloc_settings)

### Arguments

 classAlloc_settings List of inputs of the MNL model. It should contain the following. utilities: Named list of deterministic utilities . Utilities of the classes in class allocation model. Names of elements must match those in avail, if provided. avail: Named list of numeric vectors or scalars. Availabilities of classes, one element per class Names of elements must match those in classes. Values can be 0 or 1. These can be scalars or vectors (of length equal to rows in the database). A user can also specify avail=1 to indicate universal availability, or omit the setting completely. rows: Boolean vector. Consideration of which rows to include. Length equal to the number of observations (nObs), with entries equal to TRUE for rows to include, and FALSE for rows to exclude. Default is "all", equivalent to rep(TRUE, nObs). componentName: Character. Name given to model component.

### Value

The returned object depends on the value of argument functionality, which it fetches from the calling stack (see apollo_validateInputs).

• "components": Same as "estimate".

• "conditionals": Same as "estimate".

• "estimate": List of vector/matrices/arrays with the allocation probabilities for each class.

• "gradient": List containing the likelihood and gradient of the model component.

• "output": Same as "estimate".

• "prediction": Same as "estimate".

• "preprocess": Returns a list with pre-processed inputs, based on classAlloc_settings.

• "raw": Same as "estimate".

• "report": Same as "estimate".

• "shares_LL": List with probabilities for each class in an equal shares setting.

• "validate": Same as "estimate", but it also runs a set of tests to validate the function inputs.

• "zero_LL": List with probabilities for each class in an equal shares setting.

[Package apollo version 0.2.8 Index]


text:
#### Vol. 8, No. 3, 2015

 Recent Issues
 The Journal About the Journal Subscriptions Editorial Board Editors’ Interests Scientific Advantages Submission Guidelines Submission Form Editorial Login Ethics Statement Contacts Author Index To Appear ISSN: 1948-206X (e-only) ISSN: 2157-5045 (print) Other MSP Journals
Dynamics of complex-valued modified KdV solitons with applications to the stability of breathers

### Miguel A. Alejo and Claudio Muñoz

Vol. 8 (2015), No. 3, 629–674
##### Abstract

We study the long-time dynamics of complex-valued modified Korteweg–de Vries (mKdV) solitons, which are distinguished because they blow up in finite time. We establish stability properties at the ${H}^{1}$ level of regularity, uniformly away from each blow-up point. These new properties are used to prove that mKdV breathers are ${H}^{1}$-stable, improving our previous result [Comm. Math. Phys. 324:1 (2013) 233–262], where we only proved ${H}^{2}$-stability. The main new ingredient of the proof is the use of a Bäcklund transformation which relates the behavior of breathers, complex-valued solitons and small real-valued solutions of the mKdV equation. We also prove that negative energy breathers are asymptotically stable. Since we do not use any method relying on the inverse scattering transform, our proof works even under ${L}^{2}\left(ℝ\right)$ perturbations, provided a corresponding local well-posedness theory is available.

However, your active subscription may be available on Project Euclid at
https://projecteuclid.org/apde

We have not been able to recognize your IP address 34.231.247.139 as that of a subscriber to this journal.
Online access to the content of recent issues is by subscription, or purchase of single articles.

or by using our contact form.


text:
## Functional Programming on .NET part 2

7 12 2007

So most who read this aught to know what a prime number is, a number whose only factors are one and itself. This post is not on primes but does use them in a simple program which uses a simple algorithm to factor any given number. Now, many of you should know the difficulty in factoring numbers and how much of current cryptographic techniques rely on how hard the problem is. The program I wrote was written as a benchmark between the languages F# and Nemerle. Primarily as an “expressioning” benchmark and then as an aside a note on speed.

The task was to write a program which takes an unbounded positive integer and outputs a list that is the decomposition of the number into its prime factors. I decided to leverage laziness to allow me to write a flexible and short implementation. As I will get to later neither of these languages support laziness in the true sense and hence corecursion (in the sense of the dual of recursion) was unavailable, disallowing techniques such as a lazy sieve. Firstly the theory.

Theory

Given some number n, there exists primes $(p_1)(p_2)...(p_k)$ such that $n = (p_1)(p_2)...(p_k)$ and the product $(p_1)(p_2)...(p_k)$ is unique to n. This is known as the Fundamental Theorem of Arithmetic or the unique factorization theorem. I will leave the proof as an excercise to the reader. A way to get to this is to divide n by its smallest factor > 1. Store the factor and repeat on the quotient until we get to 1. Why the smallest factor? Because the smallest factor of n > 1 must be prime. Why?

Suppose that this was false and smallest factor or least divisor of n, say a, was composite. Then we would have that there exists some b and some c such that a = b * c and that b > 1 and b < a (we are working with positive integers so if b = 1 then c = a, remember that b and c < a). So b also divides n (consider n = a * d = b * c * d). But we stated earlier that a was the smallest divisor of n so by contradiction a must be prime.

So as an example we take n = 18. we have a = 9 and d = 2. 3 divides 9 so it must divide 18 and 18 / 3 = 3 * 2. Now to continue.

Suppose we keep n = 18. We know 3 is its least divisor and its prime. so p1 = 3. We do 18/3 = 6.
Next we find least divisor of 6 which is 2. so p2 = 2. Now we do 6/2 = 3. 3 is prime so 18 = 3 * 2 * 3.

So it seems to do this we need to know how to find the least divisor for a given number. Well we know that the least divisor of a number must always be prime. So if given n, to find ld of n we start from 2, check if n mod 2 = 0 and work up till we find a prime that divides n. So for 35 we’d go through 2, 3, and hit 5 as its least divisor. A way to stop unnecessary checks is to leverage the fact that if an integer n > 1 is not divisible by any primes < = $\sqrt{n}$ then n is prime.

Let n = (p)(a). Let us suppose that p is the least divisor of n, we know it must be prime. Now suppose p > $\sqrt{n}$ then a > $\sqrt{n}$ since p <= a (remember p is the smallest divisor). Then we have that $n = p \cdot a > (\sqrt{n})(\sqrt{n}) = n$ which is a contradiction since it says n > n. So at least one of p or a must be less than $\sqrt{n}$. We know p <= a so we know that at least p < $\sqrt{n}$. This proves our proposition. Note also that since $p \le a, p^2 \le p \cdot a = n$. We use that fact to lessen our checks.

The basic idea then is to define an infinite list of all natural numbers and the primes as those numbers which fail some primality test filtered out. We define that tests as any number which is its own least divisor. I will not explain the code in depth – there is plenty of documentation on both languages – this is mostly to show how easy it is to go from reasoning to implementation.

Implementation.

Below is the F# code which implements this functionality:

open System

module Primes_f =
let (+:) a alList = LazyList.cons a alList

let Z3 : bigint LazyList = LazyList.unfold (fun x -> Some(x , x + 1I)) 3I

let rec least_denominator_byfactor l n =
match l with
LazyList.Cons(p,ps) -> if n % p = 0I then
p
elif p * p > n then
n
else
least_denominator_byfactor ps n

let rec prime = function
n when n  failwith "Not positive Int"
| n when n = 1I -> false
| n when n > 1I -> least_denominator n = n
| _ -> failwith "?"

and primes =  2I +: (LazyList.filter prime Z3)

and least_denominator n = least_denominator_byfactor primes n

let rec factors = function
n when n  failwith "Not positive Int"
| n when n = 1I -> []
| n when n > 1I -> let p = least_denominator n
p :: (factors (n / p))

| _ -> failwith "?"

To display:

    let dtm = DateTime.Now
print_any ("Time is: " ^ dtm.ToString() )

//99999988888776655423101   //let z = factors 999988710I//999998888877665542310I
let rnd = Random()

let rec Generate tx l =
match tx with
0 -> l
| _ -> let z = bigint.FromInt32 (rnd.Next(1000000,int.MaxValue))
Generate (tx - 1) ( ( factors z )::l )

let list_of_primefactors =  Generate 10000 []//factors 9876543210I  //3 mins 100,000 -  20sec 10,000 .05
let dtm2 = DateTime.Now
let d = dtm2 - dtm

let primes1 = List.filter (fun l -> if List.length l = 1 then true else false) (list_of_primefactors)

let probprimes = System.Convert.ToDouble( (List.length primes1)) / 1000.0

List.map (fun c -> print_any c ; System.Console.WriteLine ("") ) primes1
print_any probprimes
//print_any list_of_primefactors

System.Console.Write ("\n\nTime is: " ^ dtm2.ToString() ^ " Took: " ^ d.ToString() ^ "\n" )

The F# code was pretty straightforward to implement. The language remains quite functional and I had a lazy list data structure which allowed semi-corecursion via unfold. Now the Nemerle code. Because of its more object-centricism (not a critique) I had to approach how I arranged the code a bit differently.

module Number
private least_denominator_byfactor(l : Primes,  n : decimal) : Decimal
def p =  l.prime

if((p==3) && (n%2 == 0))
2m
else
if(n % p == 0)
p
else if(p * p > n)
n
else
least_denominator_byfactor(l.next,n)

public least_denominator(n : Decimal) : Decimal
least_denominator_byfactor(Primes(3), n )

public IsPrime(k: Decimal) : bool
| n when n  throw ex("Not positive")
| n when n == 1 => false
| n when n == 2 => true
| n when n % 2.0m == 0 => false
| n when n > 2 => (n == least_denominator(n))

public next_prime (v : Decimal) : Decimal
def ip = IsPrime(v + 1)
match(ip)
| true => v + 1
| false => next_prime(v + 1)

public factors(number : Decimal) : list[Decimal]
| n when n  throw ex("Not Positive")
| n when n == 1 => []
| n when n > 1 => def p = least_denominator(n) ; p :: factors(n/p)

class Primes
public prime : Decimal
public next : LazyValue [Primes]

public this (p : Decimal)
prime = p
next = lazy (Primes (Number.next_prime(p)))

To display:

module Program
Main() : void

def dtm = DateTime.Now
WriteLine($"Time now is:$dtm")
// 9998876655423101
//Test Nums: 99998877665542310
def rnd = Random()

def Generate(tx, l = [])
match(tx)
| 0 => l
| _ => Generate(tx - 1, Number.factors( rnd.Next(1000000, int.MaxValue) )::l)

def list_of_primefactors = Generate(152)//Number.factors(9766427)

def dtm2 = DateTime.Now

foreach(pfactors in list_of_primefactors)
foreach(num in pfactors)
Write($"$num | ")
WriteLine("")       */

WriteLine($"Time Now is:$dtm2, Took \$(dtm2 - dtm)")

While for example in the F# code I had only functions and general functionality I had to package the code into appropriate classes in nemerle and to have behaviour sort of baked into those data types (for example I could not unfold into a generic list of primes, I had to build a next mechanism into the lazy prime class). The nemerle implementation (algorithim not language) is quite a bit less inefficient. But this is because I did not really redress the design to go with objects. For example in F# I could afford mutual recursion since F# has light weight <fast> functions and has the caveat that the lazy list caches its data.

In Nemerle the Primality check in IsPrime requires more and more objects to be created until the test is passed (note the circularity wher ld looks for a next prime which in turn calls ld to find if its prime) so it bogs down quite quickly. Indeed for small numbers (less than 10 digits) the performance is quite close – implying that the actual languages are not far apart in speed. Nonetheless the rate of increase in time usage for Nemerle is much steeper (by a factor of ~14 before change but dropping to near F#’s after for single iterations) than for F# due to the object thrashing I was doing. For example by merely injecting the following code

private ldf(k : decimal,n : decimal) : decimal
match(k)
| m when n % m == 0 => m
| m when m * m > n => n
| _ => ldf(k + 1, n)

private ld(n : decimal) : decimal
ldf(2,n)

and replacing the appropriate function in IsPrime, factoring 99998877665542310 went from taking 25 seconds to 3 seconds. A similar injection to F# increased the speed from 2 seconds to 8. F# certainly benefits from the caching in its Lazy list implementation. Both languages are eagerly evaluated except in small localized places. In effect the evaluation is more delayed than it is lazy in the true sense that one would find in Haskell. Thus for example this:

let rec sieve function = LazyList.Cons(z, zs) -> z +: sieve (LazyList.filter (fun n -> (n % z) <> 0I) zs

would fail to terminate when given sieve Z3, while similar code in Haskell knows where to stop evaluation.

Speed

As Noted earlier the F# code was fairly fast (helped, I suspect because of it caching, thus larger numbers benefited from previous computations and in having to do less work. One can really tell over multiple iterations where past iterations aid future ones. For example it took 3 minutes to factor 100,000 numbers between, 1,000,000 and 2^32. While it takes 15 secs (down from ~ 60 secs) in Nemerle to factor 120 nums, it takes F# 3 secs. In the next post I will post on Visual Studio integration, some graphs and data and improvements in the arrangement of code and primality testing for Nemerle.

One note is that F# sometimes generates confusing error messages. This error was due to using an int in place of a bigint i.e. n % 0 vs n % 0I. Error 1 Type error in the use of the overloaded operator ‘( % )’. No operator ‘%’ (full name ‘op_Modulus’) supported by the type ‘bigint’ matches the argument types. See also C:\projects\prime\prime\prime.fs(12,23)-(12,24). See also . Note the nice type system at work =)

As well, sometimes when you dont define your pattern matching cases properly it fails to build yet gives no error. Still some ways to go before being commercial. But hope it doesnt sell out when it does.

Interestingly .NET had no built in arbitrary precision integer data type (had because they just put it in for 3.5). F# does but Nemerle had to use decimal. Note though that to change I but have to rename type annotations. The code is general. Excluding the modules and classes Nemerle code and F# are very much on par.

## Functional Programming on .NET – Part 1

5 12 2007

So it seems that functional programming is the new hype. This could be a good thing but it could also be a sad/bad thing. Looking at the last major hype – object oriented programming and seeing what was lost from Smalltalk to C++ one can only hope that although concepts are brought in section by section they are done faithfully rather than as washed down ugly caricatures.

At the moment there are only 2 real offerings for functional programming on .NET and those are F# and Nemerle. There are other languages:

Scala is capable of targeting IL but I do not know its interop capabilities since it is primarily a consumer of the JVM. There is the new initiative with IronScheme but it is too early yet. As well there is the language L# but it is meant as a lightweight embeddable lisp like scripts eater. Finally there is Hugs98.NET which does not target IL but can still consume .NET code and present Haskell code for consumption. Nonetheless the lack of built in interop and the requirement to (semi)manually wrap all .NET classes to a format understood by Haskell is cumbersome. As well due to the vastly different mindsets of both platforms trying to get them to work along side each other is un-natural – like trying to get Orcs and Elves to work together.

So for now the only real contenders are F# and Nemerle. So what is it like to use each? I will not compare since that would be futile but I will note differences and such.

Firstly I will note that these languages are constrained by the platform under the hood, they are not really functional but merely compiling to .NET objects, arrays and similar. But that is not the problem (i.e. at the end of the day haskell still compiles to a bunch of gotos and alters CPU state), functional programming is more a mentality than anything. However since the functional methodologies are at the moment still not first class paradigms, to interop properly with the platform certain departures must be made from pure functional concepts in favour of a more object friendly representation. Detracting from the advantages offered by functional languages like Haskell and ML. I also suspect that the underlying “dynamic” objects (e.g. availability of reflection <- good thing?, late binding, type casting >_<) nature of the platform makes implementing a proper [static] type system a bit more of an art than necessary. Nonetheless these are all solvable issues to be addressed by Microsoft and the language designers. But at the moment there is still a fair jump down in expressive power from abstraction when going from Haskell to either of these languages.

If you go to the main Nemerle site you it will seem to be dead. One needs to look at the mailing lists to realize that the language is still being actively developed. The only detracting issue is that the language is not the focus of the developers even though it is still moving forward (amusingly, the project seems to also be very eastern with mostly eastern europe and Japan being active users). The language is very meta-powerful, allowing one to extend is capabilities using macros. Some of the things that have been done include a join calculus for concurrency, design by contract type verification support and multi-stage programming support. As a language it differs from F# as it is a bit more object oriented in the traditional sense. This I suspect may be related to why its type inference is inferior to F#’s, requiring quite a bit more type annotations. Both F# and Nemerle are on the same level of abstraction (however that would be judged) with support for succinct lightweight syntax and whitespace based formatting. Nemerle also has support for traditional ugly curly brace syntax to help bring over C# programmers. Indeed Nemerle is capable of looking very much like C# or very much like ML to suite the programmer’s tastes, making it a good candidate to act as a bridging language to functional programming.

The languages are fairly equi-capable, with the main differences between the languages being the more object oriented mentality of Nemerle (thus distracting from the functional paradigms and making its type system arguably less elegant than F#’s) and the stronger meta system of Nemerle. Finally F# is also more mature, with it being a research focus of a group of people, and production ready – it has recently been branded as an official product. Hence you are more likely to hit bugs and stumble on language inconsistencies in Nemerle than in F#. Nonetheless Nemerle is not an inferior language. It has a slightly different mentality and is an interesting combination of ideas that whose full implementation will make it a very interesting language. I feel it to be a mix of scheme, haskell and err C#.

I have split this post to reduce its size, the next post will be on a silly little lazily evaluated program in both languages I wrote which gives the prime factorization of a number. I will give my very informal speed comparisons and how each language’s mentality affected implementation.